diff --git a/.readthedocs.yml b/.readthedocs.yml index 619212a4..c2138b49 100644 --- a/.readthedocs.yml +++ b/.readthedocs.yml @@ -5,15 +5,13 @@ # Required version: 2 -# Set the version of Python and other tools you might need build: - os: ubuntu-22.04 + os: ubuntu-24.04 tools: - python: "3.10" - # You can also specify other tool versions: - # nodejs: "16" - # rust: "1.55" - # golang: "1.17" + # Capped at 3.12: docs/requirements.txt installs rabies so that + # sphinxcontrib-programoutput can run `rabies --help` at build time, and + # rabies requires Python <3.13. + python: "3.12" # Build documentation in the docs/ directory with Sphinx sphinx: @@ -26,4 +24,4 @@ formats: # Optionally declare the Python requirements required to build your docs python: install: - - requirements: docs/requirements.txt \ No newline at end of file + - requirements: docs/requirements.txt diff --git a/README.md b/README.md index 504f11a3..c54d3cbc 100644 --- a/README.md +++ b/README.md @@ -13,28 +13,26 @@ The preprocessing workflow regroups essential fMRI preprocessing steps prior to - head motion correction - susceptibility distortion correction - resampling to native or common space -- brain parcellation - slice timing correction (optional) - despiking (optional) - visual assessment of registration for quality control ### Confound correction -Following preprocessing, a range of strategies to correct fMRI confounds (e.g. motion) can then be conducted within RABIES: -- linear detrending +Following preprocessing, RABIES provides a highly-customizable confound correction pipeline regrouping several common strategies: +- frame censoring (or scrubbing) - confound regression (with several options for nuisance regressors) - frequency filtering (highpass, lowpass, bandpass) -- frame censoring (or scrubbing) - ICA-AROMA - spatial smoothing ### Analysis -Simple resting-state connectivity analyses are made available after preprocessing and confound correction. RABIES also provides a 'data diagnosis' workflow, which generates several indices of data quality and potential confounds, and conversaly, aims to improve the correction of confounds and transparency with regards to data quality: +Simple resting-state connectivity analyses are made available after preprocessing and confound correction: - seed-based functional connectivity - whole-brain connectivity matrix - group-ICA - dual regression -- data diagnosis +Importantly, the analysis stage is also complemented by a data quality assessment report (executed with `--data_diagnosis`), which generates several indices of data quality and potential confounds, with the aim of improving preprocessing strategies, data interpretation and scientific transparency. ## Notes on software design @@ -44,7 +42,7 @@ Simple resting-state connectivity analyses are made available after preprocessin - open source code - standardized input data format with [BIDS](https://bids.neuroimaging.io/) - easily shared, automatically-generated visual outputs for quality control -- containerized distribution of the software hosted on [Docker Hub](https://hub.docker.com/r/gabdesgreg/rabies) which can be downloaded via Docker and Apptainer platforms +- containerized distribution of the software through Docker/Apptainer, and conveniently installed from our [Github](https://github.com/CoBrALab/RABIES/pkgs/container/rabies) or [Docker Hub](https://hub.docker.com/r/gabdesgreg/rabies) servers. ## Citation @@ -56,12 +54,8 @@ Simple resting-state connectivity analyses are made available after preprocessin The [RABIES license](https://github.com/CoBrALab/RABIES/blob/master/LICENSE) allows for uses in academic and educational environments only. Commercial use requires a commercial license from CoBrALab , ## Acknowledgements -This software was developped by the [CoBrALab](https://cobralab.ca/), located at the Cerebral Imaging Center of the Douglas Mental Health University Institute, Montreal, Canada, in affiliation with McGill University, Montreal, Canada. This work was supported by funding from Healthy Brains, Healthy Lives (HBHL), the Fonds de recherche du Québec - Santé (FRQS) and - Nature et technologies (FRQNT), and the Natural Sciences and Engineering Research Council (NSERC) of Canada. [fMRIPrep](https://fmriprep.org/en/stable/) was an important inspirational source for this project, in particular with regards to best practices for software reproducibility and code design using Nipype. We also thank the organizers of [BrainHack School Montreal](https://school.brainhackmtl.org/), which guided the initial steps of this project in 2018. +This software was developed by the [CoBrALab](https://cobralab.ca/), located at the Cerebral Imaging Center of the Douglas Mental Health University Institute, Montreal, Canada, in affiliation with McGill University, Montreal, Canada. This work was supported by funding from Healthy Brains, Healthy Lives (HBHL), the Fonds de recherche du Québec - Santé (FRQS) and - Nature et technologies (FRQNT), and the Natural Sciences and Engineering Research Council (NSERC) of Canada. [fMRIPrep](https://fmriprep.org/en/stable/) was an important inspirational source for this project, in particular with regards to best practices for software reproducibility and code design using Nipype. We also thank the organizers of [BrainHack School Montreal](https://school.brainhackmtl.org/), which guided the initial steps of this project in 2018. ## Ask for help If you need support in using the software or experience issues that are not documented, we'll provide support on the [Github discussion](https://github.com/CoBrALab/RABIES/discussions). - -## Contributing to RABIES - -**Read our dedicated [documentation](https://rabies.readthedocs.io/en/latest/contributing.html)** diff --git a/docs/_static/custom.css b/docs/_static/custom.css new file mode 100644 index 00000000..8f40b2f7 --- /dev/null +++ b/docs/_static/custom.css @@ -0,0 +1,92 @@ +/* --------------------------------------------------------------------------- + Fully expanded sidebar navigation + + `collapse_navigation: False` makes Sphinx emit the complete toctree into + every page, but the RTD theme still hides most of it in CSS, so the nav + only opens along the branch you are currently on. See + https://github.com/readthedocs/sphinx_rtd_theme/issues/455 + + Three theme rules do the hiding: + .wy-menu-vertical li ul { display: none } + .wy-menu-vertical .toctree-lN.current .toctree-lM>ul { display: none } + .wy-menu-vertical li.current ul { display: block } + + The second is 4 classes deep, so !important is the honest way to beat it + without duplicating that whole selector chain. + --------------------------------------------------------------------------- */ +.wy-menu-vertical li ul { + display: block !important; +} + +/* With nothing left to expand, the disclosure arrows are misleading: clicking + one would appear to do nothing. */ +.wy-menu-vertical li button.toctree-expand { + display: none !important; +} + +/* The theme paints nested entries #404040 (near-black). That only works + because it shows them exclusively on the current branch, where the + background is light. Now that every branch is open, off-branch entries sit + on the dark sidebar and that colour is illegible — so give them the same + light-on-dark treatment the top-level entries get. + + `.wy-menu-vertical li.current a` ties with this selector on specificity + (both 2 classes + 2 elements), and custom.css loads later, so this would + otherwise also repaint the current branch — where the background *is* + light. The rule below puts the dark text back there. */ +.wy-menu-vertical li.toctree-l2 a, +.wy-menu-vertical li.toctree-l3 a, +.wy-menu-vertical li.toctree-l4 a { + color: #d9d9d9; +} + +.wy-menu-vertical li.current li.toctree-l2 a, +.wy-menu-vertical li.current li.toctree-l3 a, +.wy-menu-vertical li.current li.toctree-l4 a { + color: #404040; +} + +/* Indent each level and taper the size, so the hierarchy still reads at a + glance without the disclosure arrows to mark it. */ +.wy-menu-vertical li.toctree-l2 a { + padding-left: 2.427em; +} + +.wy-menu-vertical li.toctree-l3 a { + padding-left: 3.5em; + font-size: 0.85em; +} + +.wy-menu-vertical li.toctree-l4 a { + padding-left: 4.5em; + font-size: 0.8em; +} + +/* The RTD theme caps content width at 800px, which is cramped for the + parcellation tables, --help output and the wider figures. */ +.wy-nav-content { + max-width: 1000px; +} + +/* Keep wide content (help output, tables, diagrams) inside its own scroll + region rather than forcing the page to scroll horizontally. */ +.rst-content div[class^='highlight'] pre, +.rst-content table.docutils { + overflow-x: auto; +} + +/* sphinx-design cards inherit an RTD margin that double-spaces the grids. */ +.rst-content .sd-card { + margin-bottom: 0; +} + +.rst-content .sd-card .sd-card-title { + font-weight: 600; +} + +/* Definition lists are used heavily for parameter descriptions; give the term + enough weight to scan against the body text. */ +.rst-content dl:not(.docutils) > dt, +.rst-content dl.simple > dt { + font-weight: 600; +} diff --git a/docs/_static/refs.bib b/docs/_static/refs.bib index 075bedd5..faff0569 100644 --- a/docs/_static/refs.bib +++ b/docs/_static/refs.bib @@ -969,3 +969,34 @@ @ARTICLE{Biswal1995-vh keywords = "Functional connectivity;FC and Brain Networks;RABIES documentation" } + +@ARTICLE{Dorr2008-cc, + title = "High resolution three-dimensional brain atlas using an average + magnetic resonance image of 40 adult {C57Bl}/{6J} mice", + author = "Dorr, A E and Lerch, J P and Spring, S and Kabani, N and + Henkelman, R M", + journal = "Neuroimage", + volume = 42, + number = 1, + pages = "60--69", + abstract = "Detailed anatomical atlases can provide considerable interpretive + power in studies of both human and rodent neuroanatomy. Here we + describe a three-dimensional atlas of the mouse brain, manually + segmented into 62 structures, based on an average of 32 mum + isotropic resolution T(2)-weighted, within skull images of forty + 12 week old C57Bl/6J mice, scanned on a 7 T scanner. Individual + scans were normalized, registered, and averaged into one volume. + Structures within the cerebrum, cerebellum, and brainstem were + painted on each slice of the average MR image while using + simultaneous viewing of the coronal, sagittal and horizontal + orientations. The final product, which will be freely available to + the research community, provides the most detailed MR-based, + three-dimensional neuroanatomical atlas of the whole brain yet + created. The atlas is furthermore accompanied by ancillary + detailed descriptions of boundaries for each structure and + provides high quality neuroanatomical details pertinent to MR + studies using mouse models in research.", + month = aug, + year = 2008, + language = "en" +} diff --git a/docs/analysis.md b/docs/analysis.md deleted file mode 100644 index 70b71cff..00000000 --- a/docs/analysis.md +++ /dev/null @@ -1,28 +0,0 @@ -# Connectivity Analysis - -Following the completion of the confound correction workflow, RABIES allows the estimation of resting-state connectivity using standard analyses: seed-based connectivity, whole-brain connectivity, group independent component analysis (ICA) and dual regression (DR). For each analysis (except for group-ICA), RABIES will compute individualized connectivity maps for each scan separately, which can then be exported for relevant statistical analyses (e.g. group comparison) conducted outside of RABIES. - -## Correlation-based connectivity -Correlation-based analyses rely on computing a temporal correlation between different brain regions' BOLD fluctuations to estimate their functional coupling. The assessment of connectivity using correlation requires careful a priori cleaning of confounds (see [confound correction](confound_correction) and [data quality assessment](analysis_QC_target)), as various fMRI confounds introduce spurious correlations that won't be distinguished from neural activity. -(SBC_target)= -- **Seed-based connectivity** (`--seed_list`): Seed-based connectivity is the first technique developped for the mapping of connectivity during resting state {cite}`Biswal1995-vh`. The mean timecourse is first extracted from an anatomical seed of interest, and the correlation (Pearson’s r in RABIES) between this timecourse and every other voxel is computed to obtain a correlation map, representing the ‘connectivity strength’ between the seed and every other brain regions. -- **Whole-brain connectivity** (`--FC_matrix`/`--ROI_type`): This technique is an extension of the seed-based connectivity technique to encompass every brain region. That is, using the anatomical parcellation provided along with the atlas during preprocessing, the seed timecourse for every parcel is first extracted, and then the cross-correlation (Pearson’s r) is measured between every region pair. The correlation values are then re-organized into a whole-brain matrix representing the connectivity between every corresponding region pair. - -## ICA-based connectivity -The second analysis approach available within RABIES rely on the spatial decomposition of BOLD timeseries using ICA, which models the data as a linear combination of independent sources. In contrast with correlation-based connectivity, which models a single linear relationship between regions, the ICA framework accounts for multiple potentially overlapping sources of BOLD fluctuations, which may further separate confound contributions from connectivity estimates. To obtain individualized connectivity estimates, this analysis framework consists of first deriving ICA components at the group level to define sources, and then recovering individual-specific versions of the sources with dual regression {cite}`Nickerson2017-gq`. -(ICA_target)= -- **Group ICA** (`--group_ica`): RABIES uses FSL’s MELODIC ICA algorithm {cite}`Beckmann2004-yw` to derive ICA components. For group-ICA, timeseries for all scans aligned in commonspace are concatenated to group all data, before computing the ICA decomposition, yielding -$$ -Y_{concat} = A\hat{S} -$$ -where $Y_{concat}$ are the concatenated timeseries, $\hat{S}$ are the set of spatial maps defining the independent sources, and $A$ is the mixing matrix storing timecourses associated to each component. -(DR_target)= -- **Dual regression (DR)** (`--prior_maps`/`--DR_ICA`) : DR builds on the group ICA decomposition to model scan-specific versions of the group-level components, thus allowing to estimate individualized connectivity for a given brain network first identified through group-ICA {cite}`Beckmann2009-cf,Nickerson2017-gq`. DR consists of two consecutive linear regression steps, where scan-specific timecourses are first derived for each ICA component, and then a scan-specific spatial map is obtained for each component timecourse. Using multivariate linear regression (with Ordinary least square (OLS)), component timecourses are obtained with -$${\beta}_{TC} = OLS(\hat{S},Y)$$ -describing $Y = \hat{S}{\beta}_{TC} + \epsilon$ where $Y$ are the scan timeseries, $\hat{S}$ are the ICA components and ${\beta}_{TC}$ corresponds to the estimated timecourses for each component. To accurately measure connectivity amplitude in the spatial maps derived from DR, the timecourses from the first regression step must be standardized prior to the second regression{cite}`Nickerson2017-gq`. In RABIES, timecourses are thus variance-normalized using root-mean square (RMS) -$$ -{\beta}^*_{TC} = \frac{{\beta}_{TC}}{RMS({\beta}_{TC})} -$$ -where $RMS(x) = \sqrt{\frac{1}{n}\sum_{i=1}^{n}x_i^2}$. The normalized timecourses ${\beta}^*_{TC}$ are then inputed into a second regression step to derive the spatial maps ${\beta}_{SM}$ with -$${\beta}_{SM} = OLS({\beta}^*_{TC},Y^T)$$ -where $Y = {\beta}^*_{TC}{\beta}_{SM} + \epsilon$, thus completing the linear model of the timeseries. The resulting scan-specific spatial maps ${\beta}_{SM}$ will comprise information about network amplitude and shape, which may be compared across subjects or groups with further statistical tests{cite}`Nickerson2017-gq`. diff --git a/docs/analysis_QC.md b/docs/analysis_QC.md deleted file mode 100644 index 53311c7b..00000000 --- a/docs/analysis_QC.md +++ /dev/null @@ -1,56 +0,0 @@ -# Data quality assessment - -(analysis_QC_target)= - -```{toctree} ---- -maxdepth: 3 ---- -nested_docs/scan_diagnosis.md -nested_docs/distribution_plot.md -nested_docs/group_stats.md -nested_docs/optim_CR.md -``` - -Data quality can have serious impacts on analysis outcomes, leading to false findings. Rodent imaging can suffer from spurious effects on connectivity measures if potential confounds are not well accounted for, or acquisition factors, such as anesthesia levels, can influence network activity {cite}`Desrosiers-Gregoire2024-ou,Grandjean2020-fa`. To support interpretability, troubleshooting and reproducible research, RABIES includes a set of reports for conducting data quality assessment in individual scans and conducting quality control prior to network analysis at the group level. The reports are designed most specifically to evaluate the detectability of canonical brain networks and the impact of potential confounds (motion, physiological instabilities, and more). - -This page describes how to generate the reports, our guidelines for conducting quality network of network analysis, and how to include those reports in a publication. - -## Generating the reports - -At the analysis stage of the pipeline, the `--data_diagnosis` option can be selected to generate the data quality reports. To generate the report, ICA components must also be provided with `--prior_maps` and a set of components corresponding to confounds must be selected using `--conf_prior_idx` (see further details below). Connectivity can be evaluated for both dual regression and seed-based connectivity: -* **For [dual regression](DR_target)**: dual regression is always conducted using the set components from `--prior_maps`, since certain features are derived from confound components defined in `--conf_prior_idx`. On the other hand, connectivity will be evaluated in the reports for each network included in `--bold_prior_idx`. -* **For [seed-based connectivity](SBC_target)**: reports will be generated for each seed provided to `--seed_list`. However, each seed needs to be supplemented with a reference network map (a 3D Nifti file for each seed, provided with `--seed_prior_list`) which should represent the expected connectivity for the canonical network corresponding to that seed. -
-The set of reports are generated in the `data_diagnosis_datasink/` (details [here](diagnosis_datasink_target)). The interpretation of each report is described within its dedicated documentation page, and include: - -* [Spatiotemporal diagnosis](diagnosis_target): this qualitative report generated for each scan regroups a set of temporal and spatial features allowing to characterize the specific origin of data quality issues. -* [Distribution plots](dist_plot_target): quantitative report displaying the distribution of scans along measures of: specificity of network connectivity, network amplitude (for dual regression), and confound measures. Visualizing the dataset distribution can help identify outliers. -* [Group statistics](group_stats_target): for a given network, this group-level report regroups brain maps for visualizing cross-scan variability in connectivity and the group-wise correlation between connectivity and confounds. - -### Classification of group ICA components - -Ideally, the ICA components should be derived directly from the dataset analyzed by using [group ICA](ICA_target), although a [pre-computed set](https://zenodo.org/record/5118030/files/melodic_IC.nii.gz) is available by default. Newly-generated components must be visually inspected to identify the set of components corresponding to confound sources (which is inputted with `--conf_prior_idx`). This can be done by visualizing the group_melodic.ica/melodic_IC.nii.gz file, or using the automatically-generated FSL report in group_melodic.ica/report. Similarly, components corresponding to networks of interest can be identified and inputted with `--bold_prior_idx`. - -Classifying components requires careful considerations, and we recommend a conservative inclusion (i.e. not every components need to be classified, only include components which have clear feature delineating a network or a confound). Consult {cite}`Zerbi2015-nl` or {cite}`Desrosiers-Gregoire2024-ou` for more information on classifying ICA components in rodents, or the [pre-computed set](https://zenodo.org/record/5118030/files/melodic_IC.nii.gz) can be consulted as reference (the defaults for `--bold_prior_idx` and `--conf_prior_idx` correspond to the classification of these components). - -## Guidelines for analysis quality control - -![](pics/QC_framework.png) - -Below are our recommendations for how the set of quality reports can be used identify and control for the impact of data quality issues on downstream group analysis. Although the reports may be used for a breadth of applications, these guidelines are formulated most specifically for a standard resting-state fMRI design aiming to compare network connectivity between subjects or groups. In particular, the following guidelines aim to identify features of spurious or absent connectivity, remove scans where these features are prominent to avoid false results (e.g. connectivity difference is driven by motion), and determine whether these issues may confound group statistics. - -1. Inspect the [spatiotemporal diagnosis](diagnosis_target) for each scan. Particular attention should be given to the 4 main quality markers defining [categories of scan quality](quality_marker_target), and whether features of spurious or absent connectivity are prominent. -2. If spurious or absent connectivity is prominent in a subset of scans, these scans should be detected and removed to mitigate false results. This is done by setting thresholds using `--scan_QC_thresholds` for scan-level measures of network specificity and confound correlation. These measures are documented in the [distribution plots](dist_plot_target), and the specific measures for each scan ID can be consulted in the CSV file accompanying the plot. Using this CSV file, sensible threshold values can be selected for delineating scans with spurious or absent connectivity. Additionally, for dual regression analysis, `--scan_QC_thresholds` can be used to automatically detect and remove scans which present outlier values in network amplitude, which can be an indicator of spurious connectivity {cite}`Nickerson2017-gq`. By applying `--scan_QC_thresholds`, these scans won't be included for generating the group statistical report (thus the reports must re-generated after defining `--scan_QC_thresholds`). -3. Finally, the [group statistical report](group_stats_target) can be consulted to identify the main driven of variability in connectivity across scans, and whether it relates primarily to network activity or to confounds. - -If significant issues are found from this evaluation, the design of the confound correction stage may be revisited to improve quality outcomes (see [dedicated documentation](optim_CR)). - -**Disclaimer**: Although these guidelines are meant to support identifying analysis pitfalls and improve research transparency, they are not meant to be prescriptive. The judgement of the experimenter is paramount in the adopting adequate practices (e.g. network detectability may not always be expected, if studying the impact of anesthesia or inspecting a visual network in blind subjects), and the conversation surrounding what should constitute proper standards for resting-state fMRI is evolving. - -### Reporting in a publication - -All figures from the report are generated in PNG (or SVG) format, and can be shared along a publication for data transparency. Ideally, a version of the spatiotemporal diagnosis can be shared for each scan used in deriving connectivity results, together with a group statistical report and its affiliated distribution plot for each groups/datasets if the analysis involves comparing connectivity differences across subjects and/or group. -
-The set of ICA components classified as networks and confounds should be reported appropriately (e.g. melodic_IC.nii.gz file can be shared with its associated component classification). If certain scan inclusion/exclusion criteria were selected based on the quality control guidelines described above, it is particularly important to describe the observations motivating these criteria and make the associated reports readily accessible for consultation (e.g. the set of spatiotemporal diagnosis files for scans displaying spurious/absent connectivity and motivated setting a particular QC threshold with `--scan_QC_thresholds`). If the design of confound correction was defined using these tools, this should also be appropriately reported. - diff --git a/docs/bibliography.md b/docs/bibliography.md deleted file mode 100644 index 7ace31d4..00000000 --- a/docs/bibliography.md +++ /dev/null @@ -1,5 +0,0 @@ -# Bibliography - -```{bibliography} _static/refs.bib - -``` diff --git a/docs/conf.py b/docs/conf.py index 02f664ab..d3a57243 100644 --- a/docs/conf.py +++ b/docs/conf.py @@ -1,37 +1,23 @@ # Configuration file for the Sphinx documentation builder. # -# This file only contains a selection of the most common options. For a full -# list see the documentation: +# For the full list of options see: # https://www.sphinx-doc.org/en/master/usage/configuration.html -# -- Path setup -------------------------------------------------------------- - -# If extensions (or modules to document with autodoc) are in another directory, -# add these directories to sys.path here. If the directory is relative to the -# documentation root, use os.path.abspath to make it absolute, like shown here. -# -# import os -# import sys -# sys.path.insert(0, os.path.abspath('.')) - - # -- Project information ----------------------------------------------------- -project = 'RABIES Documentation' +project = 'RABIES' copyright = '2019, CoBrALab and Gabriel Desrosiers-Gregoire and Gabriel A. Devenyi and Mallar Chakravarty' author = 'CoBrALab' # The full version, including alpha/beta/rc tags release = '0.6.1' - # -- General configuration --------------------------------------------------- -# Add any Sphinx extension module names here, as strings. They can be -# extensions coming with Sphinx (named 'sphinx.ext.*') or your custom -# ones. extensions = [ "myst_parser", + "sphinx_design", + "sphinx_copybutton", "sphinx.ext.githubpages", "sphinx_rtd_dark_mode", "sphinx.ext.autosectionlabel", @@ -43,7 +29,7 @@ # to get bibliography bibtex_bibfiles = ['_static/refs.bib'] -# Choose to generate TOOD notices or not. Defaults to False +# Choose to generate TODO notices or not. Defaults to False todo_include_todos = False # Set MyST specific extensions @@ -51,14 +37,23 @@ "tasklist", "amsmath", "dollarmath", + "colon_fence", # ::: fences, so admonitions can nest and stay readable + "deflist", + "attrs_inline", ] # enable equation rendering inline myst_dmath_double_inline = True +# Generate anchors for h1-h3 so other pages can deep-link to a section +myst_heading_anchors = 3 + # Make sure the target is unique autosectionlabel_prefix_document = True +# Number and cross-reference figures, tables and code blocks +numfig = True + # Add any paths that contain templates here, relative to this directory. templates_path = ['_templates'] @@ -67,37 +62,35 @@ # This pattern also affects html_static_path and html_extra_path. exclude_patterns = ['_build', 'Thumbs.db', '.DS_Store'] - # -- Options for HTML output ------------------------------------------------- -# The theme to use for HTML and HTML Help pages. See the documentation for -# a list of builtin themes. -# -html_theme = 'groundwork' +# sphinx_rtd_dark_mode forces this theme, it is set here so the value is not +# silently disagreeing with what actually gets built. +html_theme = 'sphinx_rtd_theme' + +html_title = f'RABIES {release}' # Add any paths that contain custom static files (such as style sheets) here, # relative to this directory. They are copied after the builtin static files, # so a file named "default.css" will overwrite the builtin "default.css". html_static_path = ['_static'] +html_css_files = ['custom.css'] -# -- Options for sphinx_rtd_dark_mode ------- +# -- Options for sphinx_rtd_dark_mode ---------------------------------------- default_dark_mode = False -# Set some RTD theme config. This includes the entire navigation structure -# into the sidebar of all pages. However, expanding the sections isn't -# provided yet on the RTD theme (see -# https://github.com/readthedocs/sphinx_rtd_theme/issues/455). +# -- Options for the RTD theme ----------------------------------------------- +# The whole navigation structure is included in the sidebar of every page. html_theme_options = { 'collapse_navigation': False, 'navigation_depth': 3, + 'sticky_navigation': True, + 'titles_only': False, } -# Set some RTD theme config. This includes the entire navigation structure -# into the sidebar of all pages. However, expanding the sections isn't -# provided yet on the RTD theme (see -# https://github.com/readthedocs/sphinx_rtd_theme/issues/455). -html_theme_options = { - 'collapse_navigation': False, - 'navigation_depth': 2, -} +# -- Options for sphinx_copybutton ------------------------------------------- +# Strip shell prompts and continuation markers when copying a command +copybutton_prompt_text = r'>>> |\.\.\. |\$ ' +copybutton_prompt_is_regexp = True +copybutton_line_continuation_character = '\\' diff --git a/docs/confound_correction.md b/docs/confound_correction.md deleted file mode 100644 index 45b758c2..00000000 --- a/docs/confound_correction.md +++ /dev/null @@ -1,52 +0,0 @@ -# Confound Correction stage - -(confound_pipeline_target)= - -![Processing Schema](pics/confound_correction.png) - - -The workflow for confound correction regroups a broad set of standard tools from the human litterature. The implementation of each step is structured to follow best practices and prevent re-introduction of confounds, as recommended in {cite}`Power2014-yf` and {cite}`Lindquist2019-lq`. Importantly, each operation is optional (except detrending), and a set of operations can be selected to design a customized workflow. Optimal correction strategy can be dataset-specific, and ideally, should be tuned to address relevant quality issues identified within the dataset (see section on [data quality assessment](analysis_QC_target)). - -1. **Frame censoring** (`--frame_censoring`): Frame censoring temporal masks are derived from FD and/or DVARS thresholds, and applied first on both BOLD timeseries before any other correction step to exclude signal spikes which may bias downstream corrections, in particular, detrending, frequency filtering and confound regression{cite}`Power2014-yf`. - * Censoring with framewise displacement (see [definition](FD_target)): Apply frame censoring based on a framewise displacement threshold. The frames that exceed the given threshold, together with 1 back and 2 forward frames will be masked out{cite}`Power2012-ji`. - * Censoring with DVARS (see [definition](DVARS_target)): The DVARS values are z-scored ($DVARS_Z = \frac{DVARS-\mu}{\sigma}$, where $\mu$ is the mean DVARS across time, and $\sigma$ the standard deviation), and frames with $|DVARS_Z|>2.5$ (i.e. above 2.5 standard deviations from the mean) are removed. Z-scoring and outlier detection is repeated within the remaining frames, iteratively, until no more outlier is detected, to obtained a final set of frames post-censoring. - * `--match_number_timepoints` : This option can be selected to constrain each scan to retain the same final number of frames, to account for downstream impacts from unequal temporal degrees of freedom (tDOF) on analysis. To do so, a pre-set final number of frames is defined with `minimum_timepoint`, and a number of extra frames remaining post-censoring (taking into account edge removal in 5) ) is randomly selected and removed from the set. - -2. **Detrending** (`--detrending_order`): Linear (or quadratic) trends are removed from timeseries. Detrended timeseries $\hat{Y}$ are obtained by performing ordinary-least square (OLS) linear regression, -$$ -\beta = OLS(X,Y) -$$ -$$ -\hat{Y} = Y - X\beta -$$ -where Y is the timeseries and the predictors are $X = [intercept, time, time^2]$ ($time^2$ is included if removing quadratic trends). - -3. **ICA-AROMA** (`--ica_aroma`): Cleaning of motion-related sources using the ICA-AROMA{cite}`Pruim2015-nm` classifier. The hard-coded human priors for anatomical masking and the linear coefficients for classification were adapted from the [original code](https://github.com/maartenmennes/ICA-AROMA) to function with rodent images. ICA-AROMA is applied prior to frequency filtering to remove further effects of motion than can result in ringing after filtering{cite}`Carp2013-uf,Pruim2015-nm`. - -4. **Frequency filtering** (`--TR`/`--highpass`/`--lowpass`/`--edge_cutoff`): - 1. Simulating censored timepoints: frequency filtering requires particular considerations when applied after frame censoring, since conventional filters cannot handle missing data (censoring results in missing timepoints). To address this issue, we implemented a method described in {cite}`Power2014-yf` allowing the simulation of data points while preserving the frequency composition of the data. This method relies on an adaptation of the Lomb-Scargle periodogram, which allows estimating the frequency composition of the timeseries despite missing data points, and from that estimation, missing timepoints can be simulated while preserving the frequency profile {cite}`Mathias2004-rt`. - 2. Butterworth filter: Following the simulation, frequency filtering (highpass and/or lowpass) is applied using a 3rd-order Butterworth filter ([scipy.signal.butter](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html)). If applying highpass, it is recommended to remove 30 seconds at each end of the timeseries using `--edge_cutoff` to account for edge artefacts following filtering{cite}`Power2014-yf`. After frequency filtering, the temporal mask from censoring is re-applied to remove simulated timepoints. - -(CR_target)= -5. **Confound regression** (`--nuisance_regressors`): For each voxel timeseries, a selected set of nuisance regressors (see [regressor options](regressor_target)) are modelled using OLS linear regression and their modelled contribution to the signal is removed. Regressed timeseries $\hat{Y}$ are obtained with -$$\beta = OLS(X,Y)$$ -$$ Y_{CR} = X\beta $$ -$$ \hat{Y} = Y - Y_{CR} $$ -where $Y$ is the timeseries, $X$ is the set of nuisance timecourses and $Y_{CR}$ is the confound timeseries predicted from the model at each voxel ($Y_{CR}$ is a time by voxel 2D matrix). - -6. **Intensity scaling** (`--image_scaling`): Voxel intensity values should be scaled to improve comparability between scans/datasets. The following options are provided: - * Grand mean (**recommended**): Timeseries are divided by the mean intensity across the brain, and then multiplied by 100 to obtain percent BOLD deviations from the mean. The mean intensity of each voxel is derived from the $\beta$ coefficient from the intercept computed during **Detrending**. - * Voxelwise mean: Same as grand mean, but each voxel is independently scaled by its own mean signal. - * Global standard deviation: Timeseries are divided by the total standard deviation across all voxel timeseries. - * Voxelwise standardization: Each voxel is divided by its standard deviation. - * Homogenize variance voxelwise: if no scaling was already applied voxelwise (voxelwise mean or standardization), by selecting the option `--scale_variance_voxelwise`, timeseries are first scaled voxelwise by their standard deviation (yielding homogeneous variance distribution across voxels), and then re-scaled to preserve the original total standard deviation of the entire 4D timeseries (i.e. the global standard deviation does not change). Inhomogeneous variability distribution can be a [confound signature](quality_marker_target), thus this option may downscale their impact. `--scale_variance_voxelwise` can be applied in combination with grand mean scaling. - -7. **Smoothing** (`--smoothing_filter`): Timeseries are spatially smoothed using a Gaussian smoothing filter ([nilearn.image.smooth_img](https://nilearn.github.io/dev/modules/generated/nilearn.image.smooth_img.html)). - - -## rabies.confound_correction_pkg.confound_correction.init_confound_correction_wf [[source code](https://github.com/CoBrALab/RABIES/blob/master/rabies/confound_correction_pkg/confound_correction.py)] - -```{literalinclude} ../rabies/confound_correction_pkg/confound_correction.py -:start-after: confound_wf_head_start -:end-before: confound_wf_head_end -``` \ No newline at end of file diff --git a/docs/contributing.md b/docs/contributing.md deleted file mode 100644 index 93765cdf..00000000 --- a/docs/contributing.md +++ /dev/null @@ -1,154 +0,0 @@ -# Contributing to RABIES - -RABIES aims to provide an accessible tool responding to growing needs across the preclinical fMRI community. This effort should be community-driven, and community involvement will be paramount in achieving this goal in several respects: -- Adapting and maintaining **accessibility** for users across the broader community -- **Reproducibility and transparency**, as well as scientific scrutiny and rigor -- Defining and incorporating **best practices** across the different aspects of image processing and analysis, as well as quality control -- Leveraging appropriate expertise for the **integration of new tools** - -Suggestions for improvements can be shared using the Github [issues system](https://github.com/CoBrALab/RABIES/issues) and [discussion board](https://github.com/CoBrALab/RABIES/discussions). Contributions from developers are welcomed and encouraged. This page provides preliminary guidelines for getting started as a RABIES developer, and covers: setting a developer environment, submitting a pull request, testing and debugging, and basic instructions for adding a new module to the pipeline. We recommend discussing your proposed updates on the Github discussion board or issues prior to creating a pull request. Thank you for your support! - -## Dev environment - -For development, it is recommend to install RABIES locally, as this will make the testing and debugging process smoother. This requires installing the dependencies listed in dependencies.txt, and then installing RABIES in an appropriate python environment (e.g. using anaconda) from the Github repository. This can be done by cloning the repository, and then running ```python setup.py install```. - -### ...using a container - -It is possible to run operations using a container to avoid installing dependencies manually (however, it won't be possible to use an interface for debugging (e.g. Spyder)). This can be with `docker exec`. First, an instance of the container must be opened, which can be done by including `-d --entrypoint sh --name mycontainer` when calling `docker run`. The paths to access from the container must be set with `-v`. Here's an example: -```sh -docker run -it -v $PWD:/work_dir -v /path_to_local_RABIES_package:/RABIES:ro \ ---rm --entrypoint sh -d --name mycontainer rabies:local_testing -``` -You can then execute commands from inside the container as follows (`mycontainer` corresponds to the name set above): -```sh -docker exec mycontainer micromamba run $COMMAND -``` -To test for error, `$COMMAND` can correspond to `error_check_rabies.py --complete`. - -**Upgrading the RABIES package**: to test your updates, you must re-install RABIES inside the container. Below is a method to do so: -```sh -mkdir -p /tmp/RABIES -# copy all files from your upgraded package inside the container (/RABIES must be related to your local package with -v) -rsync -avz /RABIES/* /tmp/RABIES/. -# re-install package -cd /tmp/RABIES -python setup.py install -``` -These can be compiled into a .sh script to execute in place of `$COMMAND` above. - -## Instructions to create a pull request - -1. On github, fork the RABIES repository to have your own copy. -2. Clone your repository to carry out local modifications and testing. Use the `--recursive` option to download the submodules together with the main RABIES package. -3. Create and checkout into a new branch with `git checkout -b my_new_branch` (provide a sensible name for the branch). You are ready to make modifications to the code. -4. Testing and debugging: install your updated version of the package with ```python setup.py install```, using a proper dev environment (see above). Your can test the workflow with specific parameters by editing the ```debug_workflow.py``` script, and executing in debug mode with Spyder (see below). Before commiting changes, make sure that running ```error_check_rabies.py --complete``` completes with no error. -5. Commit and push your modifications to Github, and create a pull request from your forked repo to the original. - -## Interactive debugging with Spyder and debug_workflow.py - -Here are some recommendations for debugging using Spyder: -1. open the debug_workflow.py file in Spyder -2. find the scrips with your local installation to add breakpoints for debugging. Using `import rabies; os.path.abspath(rabies.__file__)` will provide the path to the __init__.py file of your installed package, and from there you can find file of interest and add a breakpoint where desired. -3. execute debug_workflow.py in debug mode, and run until it finds the breakpoint, and debug from there. - -## Creation of a new module and integration within a Nipype workflow - -RABIES' workflow is structured using Nipype (for more info on Nipype, see online [documentation](https://nipype.readthedocs.io/en/latest/) and [tutorial](https://miykael.github.io/nipype_tutorial/)). Preferably, a new function should be created as a Nipype interface, which has the following syntax: - -```python -from nipype.interfaces.base import ( - traits, TraitedSpec, BaseInterfaceInputSpec, - File, BaseInterface -) -class NewInterfaceInputSpec(BaseInterfaceInputSpec): - # you must select an appropriate input type with traits.type (can be Dict, File, Int, ...) - input_str = traits.Str(exists=True, mandatory=True, - desc="An input string.") - -class NewInterfaceOutputSpec(TraitedSpec): - out_file = File( - exists=True, desc="An output file.") - - -class NewInterface(BaseInterface): - """ - Describe your module. - """ - - input_spec = NewInterfaceInputSpec - output_spec = NewInterfaceOutputSpec - - def _run_interface(self, runtime): - input_str = self.inputs.input_str - - ''' - YOUR CODE - ''' - - setattr(self, 'out_file', out_file) - - return runtime - - def _list_outputs(self): - return {'out_file': getattr(self, 'out_file')} - - -``` - -You can then create a Nipype node for your interface: -```python -from .other_script import NewInterface # import your interface if from a different script -from nipype.pipeline import engine as pe - -new_interface_node = pe.Node(NewInterface(), - name='new_interface') -``` - -Instead of an interface, it is also possible to create a Nipype node from any python function: -```python -from nipype.pipeline import engine as pe -from nipype.interfaces.utility import Function - -new_function_node = pe.Node(Function(input_names=['input_1', 'input_2', ...], - output_names=['output_1', 'output_2', ...], - function=NewFunction), - name='new_function') -``` - -After creating a node which can carry the desired operation, it must be integrated within a workflow by linking up the inputs and outputs with other nodes. Below is an example of a simple workflow which conducts slice-timing correction: - -```python -from nipype.pipeline import engine as pe -from nipype.interfaces.utility import Function -from nipype.interfaces import utility as niu - -# this function creates and return a Nipype workflow which conducts slice timing correction -def init_bold_stc_wf(name='bold_stc_wf'): - - workflow = pe.Workflow(name=name) # creating a new Nipype workflow - # creating an intermediate node for storing inputs to the workflow - inputnode = pe.Node(niu.IdentityInterface( - fields=['bold_file']), name='inputnode') - # creating an intermediate node for storing outputs to the workflow - outputnode = pe.Node(niu.IdentityInterface( - fields=['stc_file']), name='outputnode') - - # preparing the node conducting STC - slice_timing_correction_node = pe.Node(Function(input_names=['in_file', 'tr', 'tpattern', 'stc_axis', - 'interp_method', 'rabies_data_type'], - output_names=[ - 'out_file'], - function=slice_timing_correction), - name='slice_timing_correction', mem_gb=1.5*opts.scale_min_memory) - - # linking up the inputnode to provide inputs to the STC node, and outputs from STC to the outputnode of the workflow - workflow.connect([ - (inputnode, slice_timing_correction_node, [('bold_file', 'in_file')]), - (slice_timing_correction_node, - outputnode, [('out_file', 'stc_file')]), - ]) - return workflow - -``` - -This example demonstrates the basic syntax of a Nipype workflow. Most likely, a new interface will be integrated as part of a pre-existing workflow (instead of creating a new one), in which case the right nodes must be linked up with the new interface. diff --git a/docs/explanation/analysis.md b/docs/explanation/analysis.md new file mode 100644 index 00000000..590eaa02 --- /dev/null +++ b/docs/explanation/analysis.md @@ -0,0 +1,124 @@ +# Connectivity analysis + +Once confound correction is complete, RABIES estimates resting-state +connectivity using standard analyses: seed-based connectivity, whole-brain +connectivity, group independent component analysis (ICA) and dual regression +(DR). + +For every analysis except group ICA, RABIES computes individualised +connectivity maps for each scan separately. These can be exported for +statistical analysis — group comparison and so on — conducted outside RABIES. + +## Correlation-based connectivity + +Correlation-based analyses compute a temporal correlation between different +brain regions' BOLD fluctuations to estimate their functional coupling. + +(SBC_target)= + +### Seed-based connectivity + +`--seed_list` + +Seed-based connectivity was the first technique developed for mapping +connectivity during rest {cite}`Biswal1995-vh`. The mean timecourse is +extracted from an anatomical seed of interest, and the correlation — Pearson's +r in RABIES — between that timecourse and every other voxel is computed, +producing a correlation map representing the connectivity strength between the +seed and every other brain region. + +### Whole-brain connectivity + +`--FC_matrix` / `--ROI_type` / `--ROI_labels_file` + +An extension of seed-based connectivity to every brain region. Using the +anatomical parcellation provided with `--ROI_labels_file`, the seed +timecourse for every parcel is extracted, then the cross-correlation (Pearson's +r) is measured between every region pair. The correlation values are +reorganised into a whole-brain matrix representing the connectivity between +every corresponding region pair. + +## ICA-based connectivity + +The second approach relies on the spatial decomposition of BOLD timeseries +using ICA, which models the data as a linear combination of independent +sources. + +Where correlation-based connectivity models a single linear relationship +between regions, the ICA framework accounts for multiple, potentially +overlapping, sources of BOLD fluctuation. This can further separate confound +contributions from connectivity estimates. + +To obtain individualised connectivity estimates, this framework first derives +ICA components at the group level to define the sources, then recovers +individual-specific versions of those sources with dual regression +{cite}`Nickerson2017-gq`. + +(ICA_target)= + +### Group ICA + +`--group_ica` + +RABIES uses FSL's MELODIC ICA algorithm {cite}`Beckmann2004-yw` to derive ICA +components. For group ICA, timeseries for all scans aligned in commonspace are +concatenated to group all data before computing the decomposition, yielding + +$$ +Y_{concat} = A\hat{S} +$$ + +where $Y_{concat}$ are the concatenated timeseries, $\hat{S}$ are the set of +spatial maps defining the independent sources, and $A$ is the mixing matrix +storing the timecourses associated with each component. + +(DR_target)= + +### Dual regression + +`--prior_maps` / `--DR_ICA` + +Dual regression builds on the group ICA decomposition to model scan-specific +versions of the group-level components, allowing individualised connectivity to +be estimated for a brain network first identified through group ICA +{cite}`Beckmann2009-cf,Nickerson2017-gq`. + +It consists of two consecutive linear regression steps. First, scan-specific +timecourses are derived for each ICA component; second, a scan-specific spatial +map is obtained for each component timecourse. + +Using multivariate OLS linear regression, component timecourses are obtained +with + +$${\beta}_{TC} = OLS(\hat{S},Y)$$ + +describing $Y = \hat{S}{\beta}_{TC} + \epsilon$, where $Y$ are the scan +timeseries, $\hat{S}$ are the ICA components and ${\beta}_{TC}$ are the +estimated timecourses for each component. + +To measure connectivity amplitude accurately in the spatial maps derived from +dual regression, the timecourses from the first regression step must be +standardised before the second regression {cite}`Nickerson2017-gq`. RABIES +variance-normalises them using root-mean square (RMS): + +$$ +{\beta}^*_{TC} = \frac{{\beta}_{TC}}{RMS({\beta}_{TC})} +$$ + +where $RMS(x) = \sqrt{\frac{1}{n}\sum_{i=1}^{n}x_i^2}$. The normalised +timecourses ${\beta}^*_{TC}$ are then fed into a second regression step to +derive the spatial maps ${\beta}_{SM}$: + +$${\beta}_{SM} = OLS({\beta}^*_{TC},Y^T)$$ + +where $Y = {\beta}^*_{TC}{\beta}_{SM} + \epsilon$, completing the linear model +of the timeseries. The resulting scan-specific spatial maps ${\beta}_{SM}$ +carry information about both network amplitude and network shape, which can be +compared across subjects or groups with further statistical tests +{cite}`Nickerson2017-gq`. + +```{seealso} +- [Analysis outputs](../reference/outputs.md#analysis-outputs) — where each result is written +- [Metric definitions](../reference/metrics.md) — precise definitions of the derived quantities +- [How to assess data quality](../how_to/assess_data_quality.md) — checking these estimates are trustworthy +``` diff --git a/docs/explanation/confound_correction.md b/docs/explanation/confound_correction.md new file mode 100644 index 00000000..0cafcf87 --- /dev/null +++ b/docs/explanation/confound_correction.md @@ -0,0 +1,187 @@ +(confound_pipeline_target)= + +# The confound correction workflow + +```{figure} ../pics/confound_correction.png +:alt: Diagram of the RABIES confound correction workflow + +The confound correction workflow. +``` + +The confound correction workflow brings together a broad set of standard tools +from the human literature. Each step's implementation follows best practices +and is structured to prevent the re-introduction of confounds, as recommended +in {cite}`Power2014-yf` and {cite}`Lindquist2019-lq`. + +Every operation is optional (at minimum the temporal mean is removed), and a set of operations can be +selected to design a customised workflow. + +```{important} +There is no universally optimal correction strategy. We provide guidelines for +tuning the pipeline to address quality issues you can +actually identify in your data — see +[Data quality assessment](data_quality.md) and +[How to optimise your confound correction strategy](../how_to/optimise_confound_correction.md). +``` + +## 1. Frame censoring + +`--frame_censoring` + +Frame censoring temporal masks are derived from FD and/or DVARS thresholds and +applied to the BOLD timeseries **first**, before any other correction step, to +exclude signal spikes which would otherwise bias detrending, frequency +filtering and confound regression {cite}`Power2014-yf`. + +Censoring with framewise displacement +: Applies frame censoring based on a [framewise displacement](FD_target) + threshold. Frames exceeding the threshold, together with 1 frame back and 2 + frames forward, are masked out {cite}`Power2012-ji`. + +Censoring with DVARS +: The [DVARS](DVARS_target) values are z-scored + ($DVARS_Z = \frac{DVARS-\mu}{\sigma}$, where $\mu$ is the mean DVARS across + time and $\sigma$ the standard deviation), and frames with $|DVARS_Z|>2.5$ + are removed. Z-scoring and outlier detection are repeated within the + remaining frames, iteratively, until no further outlier is detected. + +`--match_number_timepoints` +: Constrains every scan to retain the same final number of frames, to avoid + downstream effects of unequal temporal degrees of freedom (tDOF) on analysis. + A pre-set final number of frames is defined with `minimum_timepoint`, and the + surplus frames remaining after censoring — accounting for the edge removal in + step 4 — are selected at random and removed. + +## 2. Detrending + +`--detrending` + +Detrending is applied at the inputted polynomial order (e.g. 0 only removes the intercept, +1 for linear, 2 for quadratic, etc.). +Detrended timeseries $\hat{Y}$ are obtained by ordinary least squares (OLS) +linear regression: + +$$ +\beta = OLS(X,Y) +$$ + +$$ +\hat{Y} = Y - X\beta +$$ + +where $Y$ is the timeseries and the regressors are the polynomial expansions +of the time axis, e.g. $X = [intercept, time, time^2]$ for `--detrending order=2`. + +## 3. ICA-AROMA + +`--ica_aroma` + +Cleaning of motion-related sources using the ICA-AROMA {cite}`Pruim2015-nm` +classifier. The hard-coded human priors for anatomical masking and the linear +coefficients for classification were adapted from the +[original code](https://github.com/maartenmennes/ICA-AROMA) to function with +rodent images. + +ICA-AROMA is applied *before* frequency filtering, to remove effects of motion +that would otherwise produce ringing after filtering +{cite}`Carp2013-uf,Pruim2015-nm`. + +## 4. Frequency filtering + +`--TR` / `--highpass` / `--lowpass` / `--edge_cutoff` + +Spectral interpolation of censored timepoints +: Frequency filtering needs special handling after frame censoring, because + conventional filters cannot handle missing data. RABIES implements the method + of {cite}`Power2014-yf`, which interpolates the data while preserving its + frequency composition. It relies on an adaptation of the + Lomb-Scargle periodogram, which estimates the frequency composition of the + timeseries despite missing data points; from that estimate, missing + timepoints are simulated with the frequency profile preserved + {cite}`Mathias2004-rt`. + +Butterworth filter +: Following the simulation, highpass and/or lowpass filtering is applied using + a 3rd-order Butterworth filter + ([`scipy.signal.butter`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.butter.html)). + After filtering, the temporal mask from censoring is re-applied to remove the + interpolated timepoints. + + ```{tip} + Edge artefacts are introduced by standard frequency filters: + for a highpass at 0.01Hz, we recommend removing 30 seconds at each end of the + timeseries with `--edge_cutoff` {cite}`Power2014-yf`. + ``` + +(CR_target)= + +## 5. Nuisance regression + +`--nuisance_regressors` + +For each voxel timeseries, a selected set of +[nuisance regressors](regressor_target) is modelled using OLS linear regression +and their modelled contribution to the signal is removed. +Prior to carrying out the linear regression, a critical implementation strategy in +RABIES is to apply the same censoring, detrending and frequency filtering carried +in steps 1, 2 and 4 onto the regressors themselves to mitigate the re-introduction +of previously corrected confounds, as recommended in {cite}`Power2014-yf` and +{cite}`Lindquist2019-lq`. +After doing so, the regressed timeseries $\hat{Y}$ are obtained with + +$$\beta = OLS(X,Y)$$ + +$$ Y_{CR} = X\beta $$ + +$$ \hat{Y} = Y - Y_{CR} $$ + +where $Y$ is the timeseries, $X$ is the set of nuisance timecourses (censored, detrended +and filtered), and $Y_{CR}$ is the confound timeseries predicted from the model at each +voxel — a time-by-voxel 2D matrix. + +## 6. Intensity scaling + +`--image_scaling` + +Voxel intensity values should be scaled to improve comparability between scans +and datasets. The available options: + +Grand mean +: **Default.** Timeseries are divided by the mean intensity across the + brain, then multiplied by 100 to obtain percent BOLD deviations from the + mean. The mean intensity of each voxel is derived from the $\beta$ + coefficient of the intercept computed during **Detrending**. + +Voxelwise mean +: As grand mean, but each voxel is independently scaled by its own intercept from detrending. + +Global standard deviation +: Timeseries are divided by the total standard deviation across all voxel + timeseries. + +Voxelwise standardization +: Each voxel is divided by its own standard deviation to derive z-scored timeseries +(i.e. 0-mean and unit standard deviation). + +Homogenize variance voxelwise +: With `--scale_variance_voxelwise`, and only if no voxelwise scaling was + already applied, timeseries are first scaled voxelwise by their standard + deviation — yielding a homogeneous variance distribution across voxels — and + then re-scaled to preserve the original total standard deviation of the + entire 4D timeseries, so the global standard deviation does not change. + Inhomogeneous variability distribution can be a + [confound signature](quality_marker_target), so this option may downscale its + impact. It can be combined with grand mean scaling. + +## 7. Smoothing + +`--smoothing_filter` + +Timeseries are spatially smoothed using a Gaussian smoothing filter +([`nilearn.image.smooth_img`](https://nilearn.github.io/dev/modules/generated/nilearn.image.smooth_img.html)). + +```{seealso} +- [Nuisance regressor definitions](regressor_target) — what each regressor contains +- [`init_confound_correction_wf`](wf_confound_correction) — the source docstring +- [Confound correction outputs](../reference/outputs.md#confound-correction-outputs) +``` diff --git a/docs/explanation/data_quality.md b/docs/explanation/data_quality.md new file mode 100644 index 00000000..fc590f10 --- /dev/null +++ b/docs/explanation/data_quality.md @@ -0,0 +1,108 @@ +(analysis_QC_target)= + +# Overview of data quality assessment + +Data quality can have serious impacts on analysis outcomes, leading to false +findings. Rodent imaging can suffer from spurious effects on connectivity +measures if potential confounds are not well accounted for, and acquisition +factors such as anaesthesia level can themselves influence network activity +{cite}`Desrosiers-Gregoire2024-ou,Grandjean2020-fa`. + +To support interpretability, troubleshooting and reproducible research, RABIES +includes a set of reports for assessing data quality in individual scans and +for conducting quality control before network analysis at the group level. The +reports are designed to evaluate two main aspects: whether canonical +brain networks are detectable, and how far potential confounds — motion, +physiological instabilities, and others — have influenced the result. + +```{admonition} Where the practical instructions live +:class: seealso + +This page explains what the reports are and how they relate to one another. +For how to generate them, set inclusion thresholds and report your quality +control in a publication, see +[How to assess data quality](../how_to/assess_data_quality.md). +``` + +## The three reports + +The reports are generated by `--data_diagnosis` at the analysis stage, into +[`data_diagnosis_datasink/`](diagnosis_datasink_target). Each provides a +complementary review of the data. + +::::{grid} 1 1 3 3 +:gutter: 2 + +:::{grid-item-card} Spatiotemporal diagnosis +:link: diagnosis_target +:link-type: ref + +Qualitative, per scan. Regroups temporal and spatial features that characterise +the specific origin of a quality issue. +::: + +:::{grid-item-card} Distribution plots +:link: dist_plot_target +:link-type: ref + +Quantitative, across the dataset. Shows where each scan falls on measures of +network specificity, network amplitude and confounds — which is how outliers +become visible. +::: + +:::{grid-item-card} Group statistics +:link: group_stats_target +:link-type: ref + +Group level, per network. Brain maps of cross-scan variability in connectivity, +and the group-wise correlation between connectivity and confounds. +::: + +:::: + +## Why three levels + +```{figure} ../pics/QC_framework.png +:alt: The RABIES quality control framework + +The quality control framework. Each level conditions the validity of the next. +``` + +The three reports are not alternatives; they answer different questions, and +they depend on each other in one direction. + +The **spatiotemporal diagnosis** can help you flag most specifically what is wrong with an individual +scan — whether the signal variability carries an anatomical confound signature, +whether the network is present at all, whether network and confound timecourses +move together. It identifies the *type* of problem, which can make a +targeted correction possible. + +The **QC-FC distribution plots** turn those qualitative judgements into numbers you +can survey across the whole dataset. Doing so enables outlier detection, setting +exclusion thresholds, and detecting subtle but systematic relationships between +confounds and network measures across samples (i.e. QC-FC relationships). + +The **group statistical report** asks the question that actually matters for +common group statistical designs: is the *variability* in connectivity across scans driven by +network activity, or by confounds? Scan-level features being acceptable does +not by itself guarantee this. + +That last report depends on the first two. Either an absence of network +activity or spurious effects in a subset of scans can drive apparent network +variability, because there will be differences in the presence versus absence +of the network across scans — differences actually driven by data quality +divergences rather than by biology. This is why scan-level assumptions have to +be met before the group-level report means anything. + +## A note on judgement + +These reports and the guidelines built around them aim to identify analysis +pitfalls and improve research transparency. They are not meant to be +prescriptive. + +The judgement of the experimenter is paramount in adopting adequate practices. +Network detectability is not always expected — not when studying the impact of +anaesthesia, nor when inspecting a visual network in blind subjects. The +conversation about what should constitute proper standards for resting-state +fMRI is still evolving, and these tools are a contribution to it rather than a +settlement of it. diff --git a/docs/explanation/distribution_plot.md b/docs/explanation/distribution_plot.md new file mode 100644 index 00000000..b2c81134 --- /dev/null +++ b/docs/explanation/distribution_plot.md @@ -0,0 +1,68 @@ +(dist_plot_target)= + +# QC-FC distribution + +```{figure} ../pics/distribution_plot.png +:alt: Scatter plots of network connectivity measures against confound measures across scans + +Each point is a scan. Measures of network connectivity — specificity and +amplitude — are contrasted with measures of confounds across the sample. +``` + +The QC-FC distribution plot visualises the joint distributions of network and +confound data quality measures across the dataset, extending the per-scan qualitative +judgements from the [spatiotemporal diagnosis](diagnosis_target) into a quantitative +comparison between subjects. + +Reading the plot: + +- Points labelled in **grey** were removed using `--scan_QC_thresholds`. The + grey dotted lines are the QC thresholds selected for network specificity + (Dice overlap) and DR confound correlation. +- Among the remaining samples, and for each metric separately, scans presenting + outlier values are labelled in **orange**. Outliers are detected with a + modified Z-score threshold, set by `--outlier_threshold` and 3.5 by default. + +The derivation of each quality metric is described in the +[metric definitions](dist_plot_metrics). + +## What the report is for + +**Identify systematic QC-FC associations at the dataset-level.** Visualise the association +between network (specificity and amplitude) and a set of scan-level summary confound measures (the columns +in the plot). This complements the [group statistical report](group_stats_target) by +indicating whether a group-wise correlation in the report is +driven by a small number of outliers rather than by a dataset-wide effect. + +**Setting scan inclusion criteria.** Inspect that network specificity is +sufficient and that the temporal correlation with confounds (DR confound corr.) +is minimal, then set thresholds for scan inclusion with `--scan_QC_thresholds`. +This is the top right subplot, discussed below. + +```{seealso} +[How to assess data quality](../how_to/assess_data_quality.md) gives the +`--scan_QC_thresholds` syntax and the procedure for choosing values. +``` + + +## Inclusion criteria exemplified + +```{figure} ../pics/scan_QC_thresholds.png +:alt: Scan quality categories separated along network specificity and confound correlation axes + +Reproduced from {cite}`Desrosiers-Gregoire2024-ou`: how the +[categories of scan quality outcome](quality_marker_target) separate along +these two measures. +``` + +The measures of network specificity (Dice overlap) and temporal correlation +with confounds — where confound timecourses are extracted using the confound +components specified with `--prior_confound_idx` and measured through dual +regression — were defined in {cite}`Desrosiers-Gregoire2024-ou` for conducting +scan-level QC. + +They were selected as the measures best suited to quantifying network +detectability and spurious connectivity, and to applying inclusion thresholds +that select scans respecting the assumptions of network detectability and +minimal confound effects. + diff --git a/docs/explanation/group_statistics.md b/docs/explanation/group_statistics.md new file mode 100644 index 00000000..b2f51ad2 --- /dev/null +++ b/docs/explanation/group_statistics.md @@ -0,0 +1,72 @@ +(group_stats_target)= + +# The group statistical report + +```{figure} ../pics/group_stats_QC.svg +:alt: Group-level maps of connectivity variability and its correlation with confound measures +:width: 100% + +Group-level features of connectivity variability, for the mouse somatomotor +network. +``` + +Inspecting scan-level features is not sufficient to conclude that inter-scan +*variability* in connectivity is itself unaffected — there can be subtle +but systematic artefactual effects that impact that variability without +being easily detected from the dataset average or individual maps. +This variability is the primary driver of results in conventional +group statistical designs (e.g. comparing two different experimental groups), +in which case it is important to also assess these additional aspects of +data quality. + +## Specificity of network variability + +The standard deviation in connectivity across scans is computed voxelwise, +which visualises the spatial contrast of network variability. + +If that variability is primarily driven by network connectivity, the contrast +should reflect the anatomical extent of the network of interest, as in the +example above for the mouse somatomotor network. Otherwise it may display +spurious or absent features. For more on the development of this metric, +consult {cite}`Desrosiers-Gregoire2024-ou`. + +```{note} +**The contrast depends on sample size.** {cite}`Desrosiers-Gregoire2024-ou` +demonstrate this directly. If network connectivity is observed in individual +scans but not in this statistical report, increasing the sample size may +improve the contrast. +``` + +## Correlation with confounds + +Connectivity is correlated across subjects, at each voxel, with each of the +confound measures listed in the [metric definitions](group_QC_metrics). +This establishes how strongly connectivity is associated with potential +confounds. What constitutes a *concerning* correlation depends on the study and +on the effect size of interest: the question to ask is whether the effect size +you are looking for is much larger than the effect size of the confounds, or +comparable to it. + +## The quantitative CSV report + +A CSV file is generated alongside the figure, recording a quantitative +assessment of both aspects. The overlap between the network variability map and +the reference network map is measured using Dice overlap; for the confound +measures, the mean correlation is measured within the area of the network. See +the [group QC metric definitions](group_QC_metrics). + +```{important} +The validity of this report depends on whether the +[scan-level assumptions](dist_plot_target) of network detectability and +minimal confound effects are met. + +Either a lack of network activity or spurious effects in a subset of scans can +drive *apparent* network variability, because there will be differences in the +presence versus absence of the network across scans — but those differences +would be driven by data quality divergences rather than by biology. +``` + +```{seealso} +[How to assess data quality](../how_to/assess_data_quality.md) — the full +quality control workflow this report sits at the end of. +``` diff --git a/docs/explanation/index.md b/docs/explanation/index.md new file mode 100644 index 00000000..8121345d --- /dev/null +++ b/docs/explanation/index.md @@ -0,0 +1,39 @@ +# Explanation + +Background and discussion of how RABIES works and why it was built this way. +These pages are for understanding, not for following along — read them away +from the keyboard. + +## The processing stages + +RABIES is structured into three stages that run in sequence, each consuming the +previous one's output. + +```{toctree} +--- +maxdepth: 1 +--- +preprocessing +confound_correction +analysis +``` + +## Assessing data quality + +RABIES generates a set of reports +for characterising data quality and its impact on connectivity estimates. + +```{toctree} +--- +maxdepth: 1 +--- +data_quality +scan_diagnosis +distribution_plot +group_statistics +``` + +```{seealso} +The full methodological account is in the RABIES publication, +{cite}`Desrosiers-Gregoire2024-ou`. +``` diff --git a/docs/explanation/preprocessing.md b/docs/explanation/preprocessing.md new file mode 100644 index 00000000..39042a2b --- /dev/null +++ b/docs/explanation/preprocessing.md @@ -0,0 +1,99 @@ +# The preprocessing workflow + +```{figure} ../pics/preprocessing.png +:alt: Diagram of the RABIES preprocessing workflow + +The RABIES preprocessing workflow. +``` + +Preprocessing fMRI scans prior to analysis requires, at minimum, anatomical +alignment of scans to a common space, head realignment to correct for motion, +and correction of the susceptibility distortions arising from the echo-planar +imaging (EPI) acquisition of functional scans. The core preprocessing pipeline +in RABIES carries out each of these with state-of-the-art processing tools and +techniques. + +## Alignment to the common space + +Structural images, acquired alongside the EPI scans, are initially corrected +for inhomogeneities (**Structural inhomogeneity correction**) and then +registered together so that different MRI acquisitions can be aligned. + +That registration works by generating an unbiased, data-driven template +(**Unbiased template generation**) through the iterative non-linear +registration of each image to the dataset consensus average, where the average +is updated at each iteration to provide an increasingly representative dataset +template ([optimized_antsMultivariateTemplateConstruction](https://github.com/CoBrALab/optimized_antsMultivariateTemplateConstruction); +{cite}`Avants2011-av`). + +The finalised template, after the last iteration, provides a representative +alignment of each MRI session to a template sharing the acquisition properties +of the dataset — brain shape, field of view, anatomical contrast — which makes +it a stable registration target for cross-subject alignment. This +newly-generated unbiased template is then itself registered to an external +reference atlas (**Atlas registration**), which supplies both an anatomical +segmentation and a common space comparable across studies. + +```{note} +This is why RABIES builds a study-specific template rather than registering +every scan directly to the atlas: the intermediate target resembles your data, +so each individual registration has less work to do and is less likely to +fail. The cost is computation time, which is why `fast_commonspace=true` +exists to skip it. +``` + +## EPI motion and distortion corrections + +A volumetric EPI image is first derived using a trimmed mean across the EPI +frames, after an initial motion realignment step (**3D EPI generation**). Using +this volumetric EPI as a target, the head motion parameters are estimated by +realigning each EPI frame to the target with a rigid registration +(**Head motion estimation**). + +To correct EPI susceptibility distortions, the volumetric EPI is first +subjected to an inhomogeneity correction step +(**Functional inhomogeneity correction**), then registered non-linearly to the +anatomical scan from the same MRI session, which yields the geometrical +transforms required to recover brain anatomy {cite}`Wang2017-ci` +(**Susceptibility distortion estimation**). + +## Derivation of preprocessed EPI timeseries + +The transforms providing head motion correction, susceptibility distortions +and alignment to the common space are concatenated into a single resampling operation +— avoiding multiple resampling — applied at each EPI frame {cite}`Esteban2019-rs` (**Frame-wise resampling**). +This generates the preprocessed EPI timeseries in common space, while alternatively +the transforms to common space can be dropped to generate instead native space timeseries +using the `--resampling_space` parameter. + +```{important} +Concatenating the transforms matters. Resampling an image is lossy, so +applying motion correction and then distortion correction as two separate +resampling steps blurs the data twice. RABIES composes the transforms first +and resamples once. +``` + +## Working without structural scans + +Structural scans are recommended but not required. An alternative workflow, +selected with `--bold_only`, preprocesses an input dataset containing only EPI +functional images. + +In this workflow the volumetric EPI corrected during +**Functional inhomogeneity correction** replaces the structural image for the +purpose of common space alignment, and is used to generate the unbiased +template, which is in turn registered to the reference atlas. This final +registration to the atlas accounts for the estimation of susceptibility +distortions, in place of the registration to a same-session structural image. + +```{note} +When using the RABIES default mouse atlas, `--bold_only` also switches the +default template to an EPI reference template, which is a more robust target +for EPI registration than a structural reference. +``` + +```{seealso} +- [Preprocessing QC outputs](../reference/qc_outputs.md) — how to verify each step succeeded +- [How to troubleshoot registration](../how_to/troubleshoot_registration.md) — what to do when one did not +- [Workflow reference](../reference/workflows.md) — the source docstrings for every module above +``` diff --git a/docs/explanation/scan_diagnosis.md b/docs/explanation/scan_diagnosis.md new file mode 100644 index 00000000..6f6466f9 --- /dev/null +++ b/docs/explanation/scan_diagnosis.md @@ -0,0 +1,189 @@ +(diagnosis_target)= + +# The spatiotemporal diagnosis + +The spatiotemporal diagnosis is a visual report generated for each scan +independently, after dual regression or seed-based connectivity analysis. It +displays a large set of temporal and spatial features for that scan, supporting +assessment of potential data quality issues and of whether network connectivity +has been affected. + +It is generated by `--data_diagnosis` at the analysis stage, into +`data_diagnosis_datasink/figure_temporal_diagnosis/` and +`data_diagnosis_datasink/figure_spatial_diagnosis/`. + +```{note} +Unless stated otherwise, all metrics are computed from fMRI timeseries **after** +the confound correction stage. +``` + +## An example report + +```{figure} ../pics/spatiotemporal_diagnosis.png +:alt: The spatiotemporal diagnosis report for a scan with clear network connectivity + +A scan with little confound signature and clear network connectivity. Each +spatial map is shown along 6 cortical slices, overlapped onto the anatomical +template in common space. +``` + +In this example both dual regression and seed-based connectivity were +conducted: DR network 0 and SBC network 1 correspond to analysis of the +somatomotor network, while DR network 1 and SBC network 0 correspond to the +default mode network. Network maps from dual regression (DR) or seed-based +connectivity (SBC) are thresholded to include the top X% of voxels with the +highest values, X% being defined by `--brainmap_percent_threshold`. + +The features are detailed below. For the exact computation behind each one, see +[Metric definitions](metrics_target). + +### Temporal features, top to bottom + +**Power spectrum** +: The frequency power spectrum, showing the dominant frequency profile. Always + shown averaged across the whole brain; spectra for individual seeds can be + added with `--plot_seed_frequencies`. + +**Carpet plot** +: The entire fMRI timeseries as a time-by-voxel 2D matrix. This makes global + fluctuations in signal intensity visible, which is a proxy for various global + artefacts {cite}`Power2017-wn`. Censored frames are labelled in red. + +**Translation and rotation head motion parameters** +: The 6 rigid body parameters estimated during preprocessing, tracking head + position across the scan duration. + +**Framewise distance from the mean image** +: The mean square error between each frame and the average EPI, computed as the + tri-mean across time voxelwise, and computed **prior to** confound correction + to map changes in brain contrast over time. This complements framewise + displacement, since it captures changes in raw image contrast across time + that remain after head motion realignment and might not be well captured by + framewise displacement. The metric is also output by `--hmc_qc_report` at the + `preprocess` stage, and can be used for censoring via `--frame_censoring` at + the `confound_correction` stage. + +**Framewise displacement and DVARS** {cite}`Power2012-ji` +: Framewise displacement quantifies movement between consecutive frames, + revealing the timing and amplitude of spontaneous motion. DVARS — the + root-mean-square of the timeseries' temporal derivative — reveals shifts in + global fMRI signal intensity, which can also indicate suspicious spikes in + signal. + +**Variance captured and removed by confound regression** +: The linear model fitted during nuisance regression provides an estimate of + the nuisance timeseries, from which the variance removed at each frame can be + quantified. The report displays the absolute variance, as the RMS of the + nuisance timeseries across voxels for each frame, and the relative variance + removed, as the $R^2$ at each frame. Together these visualise confound + effects and show whether confound regression appropriately modelled the + confounds detected in the other temporal features. + +**Mean signals within anatomical masks** +: The mean timecourse within the whole-brain, edge, white matter and CSF masks. + The edge, WM and CSF regions are susceptible to motion and physiological + confounds {cite}`Pruim2015-nm`, and contrasting them with the whole-brain + global signal helps assess artefactual effects across the brain. + +**Mean amplitude of network vs confound timecourses** +: The averaged timecourses for network analyses and for confound sources are + compared, to assess whether network amplitude is spurious — that is, + correlated with the confound timecourse. Confound timecourses are modelled by + running dual regression with a complete set of ICA components representing a + mixture of networks and confounds of various origins, and compiling the + timecourses from the confound components. RABIES uses + [this set](https://zenodo.org/records/19069284/files/melodic_IC.nii.gz) of + ICA components for mice by default. + +### Spatial features, top to bottom + +**Spatial distribution in signal variability (BOLDSD)** +: The signal variability, i.e. standard deviation, at each voxel. This map + indexes whether significant confounds are contributing to the signal. Without + influence from confounds, as in the example above, signal variability is + largely homogeneous. + +**Confound regression variance explained (CRSD)** +: The variance explained by confound regression at each voxel, taken as the + standard deviation of the modelled confound timecourse. This contrasts the + amplitude of confound effects spatially, and can specifically delineate the + presence of confounds and identify their type. In the example, minor motion + signatures are identified. + +**Confound regression variance explained proportion** +: As CRSD, but showing the proportion of variance explained ($R^2$). + +**Global signal covariance** +: The covariance of each voxel with the global signal. The contrast in this map + reveals the predominant source of global signal fluctuation, which takes + various forms depending on the contributions from neural networks and from + confounds. In the ideal case there is predominant contrast in grey matter, + with a shape reminiscent of a brain network, as in the example above. + +**Network spatial maps** +: The spatial network maps fitted using dual regression, or seed-based + analysis, for the selected networks of interest — in the example, the + somatomotor and default mode networks. These fits give insight into the + quality of the network analysis and how it may affect downstream statistics. + +```{note} +$CR_{SD}$ and $CR_{R^2}$ are computed from the regressors specified with +`--nuisance_regressors` during confound correction. If no regressors are +specified, they are still estimated using regression on the 6 motion +parameters — but that regression is not applied to remove signal from the +timeseries. +``` + +(quality_marker_target)= + +## Four categories of quality outcome + +```{figure} ../pics/diagnosis_key_markers.png +:alt: Five key diagnostic features enable identifying four categories of scan quality + +Five key diagnostic features enable identifying four categories of scan quality. +``` + +A subset of the features above are the most decisive in determining scan +quality in relation to connectivity analysis. We exemplify here how these +can distinguish between 4 main categories of scan quality. + +**BOLD variability** +: Presents a homogeneous contrast in uncorrupted scans, and otherwise reveals + the anatomical signature of a confound — this can allow recognising the type of + confound present. + +**Global signal covariance** +: Sensitive to both non-neural confounds (the *spurious* category) and network + signatures (the *specific* category). It therefore reflects whether network + or confound sources dominate coordinated fluctuations, and delineates the + most likely contributors to downstream connectivity measures. + +**Network map** +: By visualising the network maps, we can directly intuit whether they + reproduce expected anatomical features of the network of interest. + This is what establishes that the network + is not absent (the *absent* category), and that its shape is not distorted by + spurious features (the *spurious* category). + +**Network and confound timecourses** +: The respective timecourses reveal relationships between network + amplitude and confounds in the temporal domain. This does not describe the + *type* of confound, but a strong correlation in the temporal domain + is the most conclusive indicator of spurious connectivity, + and an important complement to inspecting network shape — + spurious effects can affect amplitude while leaving shape largely intact. + +These 4 features are sufficient to capture the essential characteristics of +network detectability and spurious connectivity at the single scan level. The +remaining features provide additional detail about timeseries properties, +motion parameters and confound regression, and further support characterising +the specific origin of a confound — for example, establishing that a +correlation between network and confound timecourses originates from framewise +displacement, i.e. from motion. + +```{seealso} +- [How to assess data quality](../how_to/assess_data_quality.md) — acting on what you see here +- [Distribution plot](distribution_plot.md) — turning these judgements into thresholds +- [Metric definitions](metrics_target) — the computation behind each feature +``` diff --git a/docs/faq.md b/docs/faq.md deleted file mode 100644 index 9e4c2800..00000000 --- a/docs/faq.md +++ /dev/null @@ -1,55 +0,0 @@ -# FAQ - -## How do I run minimal preprocessing? - -Some users arrive with functional data that was already preprocessed with their own pipeline, and only want to use RABIES for confound correction and analysis. - -Running `rabies confound_correction` directly on such data fails with missing `.pkl` file errors. Those `.pkl` files are serialized states of the internal preprocessing workflow and cannot be written by hand. Even with the `.pkl` in place, the confound correction stage expects the full range of output files produced by the preprocessing stage, so missing file errors would follow. - -The practical solution is to run a *SHAM* preprocessing, in which the correction and registration steps are turned off. RABIES still computes the intermediary outputs that the confound correction and analysis stages require, while leaving the image data largely unchanged. It is not a strict pass-through, however, and the operations which remain applied are described under *What still happens* below: - -```sh -rabies preprocess bids_inputs/ preprocess_outputs/ \ - --anat_inho_cor method=disable \ - --bold_inho_cor method=disable \ - --commonspace_reg template_registration=no_reg,fast_commonspace=true \ - --bold2anat_coreg registration=no_reg \ - --no_HMC -``` - -| Parameter | What it turns off | -|---|---| -| `--anat_inho_cor method=disable` | inhomogeneity correction of the structural images | -| `--bold_inho_cor method=disable` | inhomogeneity correction of the EPI images | -| `--commonspace_reg template_registration=no_reg,fast_commonspace=true` | unbiased template generation, and the registration to the reference atlas (an identity transform is used instead) | -| `--bold2anat_coreg registration=no_reg` | EPI-to-structural coregistration, i.e. the susceptibility distortion correction | -| `--no_HMC` | the *application* of head motion correction to the resampled timeseries. Head motion parameters are still estimated and remain available to `--conf_list`, `--frame_censoring` and `--data_diagnosis` | - -If your dataset has no structural scans, add `--bold_only`, in which case `--anat_inho_cor` and `--bold2anat_coreg` no longer apply. - -### Caveat: your data must already be in commonspace - -Setting `template_registration=no_reg` does not skip the resampling to commonspace, it replaces the estimated transform with an identity transform. The commonspace outputs are therefore only meaningful if your input images already overlap with the template provided through `--anat_template`, which is the file that defines the commonspace. If they do not overlap, the commonspace timeseries and the atlas masks (`--brain_mask`, `--WM_mask`, `--CSF_mask`, `--vascular_mask`) applied downstream will not correspond to your data. - -Image orientation is worth confirming before you start, see [Troubleshooting](troubleshooting.md). - -### What still happens - -Operations which alter the data but are off by default stay off, and should not be added: `--apply_STC`, `--apply_despiking`, `--detect_dummy`, `--log_transform`, `--anat_autobox`, `--bold_autobox` and `--oblique2card`. - -The timeseries are still resampled onto the output grid, using the identity transforms described above. Use `--commonspace_resampling` and `--anatomical_resampling` to control the output voxel dimensions, and `--interpolation` to select the interpolator. - -Two further operations are applied unconditionally and cannot be turned off: - -* **Negative values are clipped to zero** when the preprocessed timeseries are written out. If your data legitimately contains negative values, for instance because it was already demeaned or detrended by your own pipeline, those voxels will be set to zero. Bring in data on a positive scale, and leave centering to `--detrending` at the confound correction stage. -* **The output is cast** to the type given by `--data_type`, which is `float32` by default. - -A SHAM preprocessing is therefore a minimal pass, not a strictly non-modifying one. - -### Alternative: `--read_datasink` - -`rabies confound_correction --read_datasink` reads the preprocessing outputs from the datasink folders rather than from the saved workflow graph, which removes the need for the `.pkl` file. This requires reproducing the RABIES preprocessing output structure and file naming exactly, and is generally more work than running a SHAM preprocessing. - -### Still not covered? - -If your use case needs settings that are not exposed, open a [discussion](https://github.com/CoBrALab/RABIES/discussions) describing explicitly what you need. Providing example data lets us work out an implementation supporting your use case. diff --git a/docs/how_to/assess_data_quality.md b/docs/how_to/assess_data_quality.md new file mode 100644 index 00000000..9a332fd0 --- /dev/null +++ b/docs/how_to/assess_data_quality.md @@ -0,0 +1,179 @@ +# How to assess data quality + +This guide covers generating the RABIES data quality reports, using them to +decide which scans to keep, and reporting what you did in a publication. For +what the reports mean and why they exist, see +[Data quality assessment](../explanation/data_quality.md). + +The guidance below is written for a standard resting-state fMRI design in which +you compare network connectivity between subjects or groups. The aim is to +identify features of spurious or absent connectivity, remove the scans where +those features dominate, and establish whether the remaining issues confound +your group statistics. + +## Generate the reports + +Pass `--data_diagnosis` at the analysis stage. It needs a set of ICA components +via `--prior_maps`, with the components corresponding to confounds identified +through `--prior_confound_idx`: + +```sh +rabies -p MultiProc analysis confound_correction_outputs/ analysis_outputs/ \ + --data_diagnosis \ + --prior_maps melodic_IC.nii.gz \ + --prior_bold_idx 5 12 19 \ + --prior_confound_idx 0 1 2 6 7 8 \ + --DR_ICA +``` + +The reports appear in +[`data_diagnosis_datasink/`](../reference/outputs.md#analysis-outputs). + +Connectivity can be evaluated for either analysis, or both: + +For [dual regression](DR_target) +: Dual regression is always run using the full set of components from + `--prior_maps`, because several report features are derived from the confound + components named in `--prior_confound_idx`. Connectivity itself is evaluated for + each network listed in `--prior_bold_idx`. + +For [seed-based connectivity](SBC_target) +: A report is generated for each seed given to `--seed_list`. Each seed must be + accompanied by a reference network map — a 3D NIfTI file per seed, passed + through `--seed_prior_list` — representing the connectivity expected for the + canonical network that seed belongs to. + +### Classify your group ICA components + +Ideally the components come from the dataset you are analysing, derived with +[group ICA](ICA_target). A +[pre-computed set](https://zenodo.org/records/19069284/files/melodic_IC.nii.gz) +for mice is used by default. + +Newly generated components must be inspected visually to identify which +correspond to confound sources. Visualise `group_melodic.ica/melodic_IC.nii.gz`, +or use the FSL report generated automatically in `group_melodic.ica/report`. +Pass the confound components to `--prior_confound_idx` and the networks of interest +to `--prior_bold_idx`. + +```{tip} +Classify conservatively. Not every component needs a label — include only +those with a clear feature delineating a network or a confound. The defaults +for `--prior_bold_idx` and `--prior_confound_idx` correspond to the classification +of the pre-computed set, which you can consult as a reference. + +For guidance on classifying ICA components in rodents, see +{cite}`Zerbi2015-nl` and {cite}`Desrosiers-Gregoire2024-ou`. +``` + +## Work through the reports + +```{figure} ../pics/QC_framework.png +:alt: The RABIES quality control framework, from scan-level diagnosis to group statistics + +The quality control framework: scan-level diagnosis feeds scan inclusion +decisions, which in turn condition the validity of the group-level report. +``` + +### 1. Inspect each scan + +Read the [spatiotemporal diagnosis](diagnosis_target) for every scan. Pay +particular attention to the four main quality markers that define the +[categories of scan quality](quality_marker_target), and judge whether features +of spurious or absent connectivity are prominent. + +### 2. Remove scans with spurious or absent connectivity + +If those features are prominent in a subset of scans, remove those scans to +mitigate false results. Set thresholds with `--scan_QC_thresholds` on the +scan-level measures of network specificity and confound correlation: + +```sh +rabies -p MultiProc analysis confound_correction_outputs/ analysis_outputs/ \ + --data_diagnosis \ + --prior_maps melodic_IC.nii.gz \ + --prior_bold_idx 5 12 19 --prior_confound_idx 0 1 2 6 7 8 --DR_ICA \ + --scan_QC_thresholds '{DR:{Dice:[0.3,0.3,0.3],Conf:[0.25,0.25,0.25],Amp:false}}' +``` + +The value is a dictionary expression, quoted so the shell leaves it alone. Per +analysis (`DR`, `SBC` or `NPR`) you can set: + +`Dice` +: Minimum network detectability, as Dice overlap with the prior. A list of + values between 0 and 1, matched in order to `--prior_bold_idx` for DR and + NPR, or to `--seed_list` for SBC. Either give an empty list, or give exactly + as many thresholds as there are networks. + +`Conf` +: Maximum temporal correlation with the dual regression confound timecourses. + Same list rules as `Dice`. + +`Amp` +: `true` to automatically remove scans with outlier network amplitude, which + can indicate spurious connectivity {cite}`Nickerson2017-gq`. + +**Sensible threshold values** should be selected by relating scans flagged in +step 1 (i.e. those that present spurious/absent features in the [spatiotemporal diagnosis](diagnosis_target)) +to their associated `Dice`/`Conf` values listed in the [distribution plots](dist_plot_target) +and the accompanying CSV file, which gives the measures per scan ID. +We do not recommend blindly applying a threshold value listed in a previous +publication, as legitimate `Dice`/`Conf` values will differ depending +on the image signal-to-noise ratio and/or preprocessing decision +(e.g. applying a lowpass filter will systematically increase `Conf` values). + +```{important} +Scans excluded by `--scan_QC_thresholds` are excluded from the group +statistical report, so it is best to regenerate these reports after +you set the thresholds. +``` + +### 3. Check the group level + +Consult the [group statistical report](group_stats_target) to identify the main +driver of connectivity variability across scans, and whether there +are systematic group-level associations with confound metrics. + +### 4. Revisit confound correction if needed + +If significant issues remain, consider redesigning the confound correction +stage — see [How to optimise your confound correction strategy](optimise_confound_correction.md). + +```{admonition} These guidelines are not prescriptive +:class: caution + +They are meant to support identifying analysis pitfalls and improving research +transparency. The judgement of the experimenter is paramount: network +detectability is not always expected, for instance when studying the impact of +anaesthesia or inspecting a visual network in blind subjects. The conversation +about what should constitute proper standards for resting-state fMRI is still +evolving. +``` + +## Report your quality control in a publication + +A central motivation for implementing these automatically-generated +quality reports is to encourage and improve scientific transparency and +study comparison. Every figure in the report is generated as PNG or SVG +and can be shared alongside a publication. + +We specifically recommend sharing: + +- The spatiotemporal diagnosis report from each scan used to derive connectivity results. +- A group statistical report and its affiliated distribution plot, for each + group or dataset, if the analysis compares connectivity across subjects + and/or groups. +- The set of ICA components classified as networks and confounds — for example + the `melodic_IC.nii.gz` file, with its component classification. + +If you excluded scans based on the guidelines above, describe the observations +that motivated your criteria and make the associated reports accessible: the +spatiotemporal diagnoses for the scans that displayed spurious or absent +connectivity and motivated a particular `--scan_QC_thresholds` value. If you +designed your confound correction using these tools, report that too. + +```{seealso} +- [Data quality assessment](../explanation/data_quality.md) — what these reports are for +- [Metric definitions](../reference/metrics.md) — how every quantity is computed +- [Analysis outputs](../reference/outputs.md#analysis-outputs) — where each file lands +``` diff --git a/docs/how_to/change_template.md b/docs/how_to/change_template.md new file mode 100644 index 00000000..0edd86d5 --- /dev/null +++ b/docs/how_to/change_template.md @@ -0,0 +1,27 @@ +# How to override the default common space template + +The files that will define the common space template and associated masks are +controlled by the following parameters at the `preprocess` stage: +- `--anat_template`: the structural image that provides the target for common space + alignment. +- `--brain_mask`: the brain mask file. +- `--WM_mask`: a mask for the white matter. +- `--CSF_mask`: a mask for the cerebrospinal fluid canals. +- `--vascular_mask`: a mask labelling major blood vessels. + +RABIES will automatically input a file for each parameter by default (see [built-in commonspace template and atlas](../reference/template_files.md)), +but you may want to choose a separate template for your study (for instance if you are not inputting adult mouse data). To do so, you will need +to provide the full path to an adequate set of NIfTI-formatted files. At minimum, both `--anat_template` and `--brain_mask` must be provided, +as they are key components of the registration pipeline. The other brain parcellations are optional, but certain downstream options for +confound correction or generating QC reports will be disabled. + +```{important} +The template file `--anat_template` defines the common space coordinates for the entire pipeline. All other +mask files are expected to overlap with this template file, as assessed from the metadata coordinates of the +NIfTI header. An error will be thrown if a mismatch is detected. +``` + +```{important} +At minimum, the image orientation of the template and input data must match. It is recommended to follow the +RAS+ orientation convention that RABIES inspects, see [How to check image orientation](./check_orientation.md). +``` diff --git a/docs/how_to/check_orientation.md b/docs/how_to/check_orientation.md new file mode 100644 index 00000000..9c82b72b --- /dev/null +++ b/docs/how_to/check_orientation.md @@ -0,0 +1,57 @@ +# How to check image orientation + +RABIES assumes input data is oriented according to the NIfTI standard (RAS+, +i.e. Right–Anterior–Superior). Incorrectly oriented images are one of the most +common causes of registration failure and of results that look wrong for no +apparent reason, because the ANTs/ITK tools RABIES calls will happily register +a mis-oriented brain to the template and produce nonsense. + +```{important} +Check orientation before you report a bug, and before you start +troubleshooting anything else. It costs five minutes and rules out the most +likely explanation. +``` + +## Check with ITK-SNAP + +[ITK-SNAP](https://www.itksnap.org/pmwiki/pmwiki.php) is a free, open-source +medical image viewer that reads NIfTI orientation information correctly. Other +viewers do not always, which is why it is the recommended tool here. + +1. **Open your image.** *File → Open Main Image…*, and load your NIfTI + anatomical scan. + +2. **Check the anatomy is where you expect it.** ITK-SNAP shows axial, coronal + and sagittal views. The nose should be anterior, the top of the head + superior. + +3. **Compare against your reference template.** Open the atlas you intend to + pass to `--anat_template` (SIGMA, DSURQE, Fischer rat, whichever applies) in + a second ITK-SNAP window. Structures should appear in similar positions, and + the orientation labels should match. + +4. **Check the orientation labels.** Each view is labelled **R/L**, **A/P** and + **S/I**. Confirm each corresponds to the real anatomical direction in your + scan, then move the cursor: the crosshair should track consistently across + all three views, so that dragging towards the label **R** moves towards + anatomical right. + +If any of these disagree, the orientation stored in your NIfTI header does not +describe your data, and you need to fix it before running RABIES. This is +usually introduced during conversion from the scanner format. + +```{seealso} +For Bruker data, [BrkRaw](https://brkraw.github.io/) handles the raw-to-NIfTI +conversion; the CoBrALab maintains +[notes on the conversion](https://github.com/CoBrALab/documentation/wiki/bruker2nifti-conversion). +``` + +## If orientation is correct and registration still fails + +Move on to [How to troubleshoot registration](troubleshoot_registration.md), +which covers tuning inhomogeneity correction, masking and registration stages. +Read the [preprocessing QC report](../reference/qc_outputs.md) first, to +establish which registration step is the one that failed. + +When you do report a problem, follow the +[issue template](https://github.com/CoBrALab/RABIES/blob/master/.github/ISSUE_TEMPLATE/standard-bug-report.md). diff --git a/docs/how_to/contribute.md b/docs/how_to/contribute.md new file mode 100644 index 00000000..d746442d --- /dev/null +++ b/docs/how_to/contribute.md @@ -0,0 +1,218 @@ +# How to contribute to RABIES + +RABIES aims to provide an accessible tool responding to growing needs across +the preclinical fMRI community. This effort should be community-driven, and +community involvement is paramount in several respects: + +- adapting and maintaining **accessibility** for users across the broader community +- **reproducibility and transparency**, as well as scientific scrutiny and rigour +- defining and incorporating **best practices** across image processing, analysis and quality control +- leveraging appropriate expertise for the **integration of new tools** + +Suggestions for improvements can be shared through the GitHub +[issues system](https://github.com/CoBrALab/RABIES/issues) and +[discussion board](https://github.com/CoBrALab/RABIES/discussions). + +## Set up a development environment + +Install RABIES locally rather than working through a container — testing and +debugging are much smoother that way. Install the dependencies listed in +[`dependencies.txt`](https://github.com/CoBrALab/RABIES/blob/master/dependencies.txt), +then install RABIES from a clone of the repository into a Python environment of +your choice: + +```sh +git clone --recursive https://github.com/CoBrALab/RABIES.git +cd RABIES +python setup.py install +``` + +```{note} +Use `--recursive`. RABIES pulls in submodules, and a clone without them will +not run. +``` + +### Working inside a container instead + +You can run operations in a container to avoid installing dependencies by hand, +at the cost of losing interactive debugging (Spyder and similar will not be +available). + +Open a persistent container instance with `-d --entrypoint sh --name mycontainer`, +binding the paths you need with `-v`: + +```sh +docker run -it -v $PWD:/work_dir -v /path_to_local_RABIES_package:/RABIES:ro \ + --rm --entrypoint sh -d --name mycontainer rabies:local_testing +``` + +Then execute commands inside it: + +```sh +docker exec mycontainer micromamba run $COMMAND +``` + +To check for errors, `$COMMAND` can be `error_check_rabies.py --complete`. + +To test your changes you must reinstall RABIES inside the container: + +```{code-block} sh +:caption: Reinstalling the package inside a running container + +mkdir -p /tmp/RABIES +# copy the upgraded package into the container +# (/RABIES must be bound to your local package with -v) +rsync -avz /RABIES/* /tmp/RABIES/. +cd /tmp/RABIES +python setup.py install +``` + +Compile these into a `.sh` script and run it in place of `$COMMAND` above. + +## Submit a pull request + +1. **Fork** the RABIES repository on GitHub. +2. **Clone** your fork, with `--recursive` to pull the submodules. +3. **Branch**: `git checkout -b my_new_branch`, with a name describing the + change. You are ready to modify the code. +4. **Test and debug.** Install your updated package with + `python setup.py install` in a proper development environment. Test the + workflow with specific parameters by editing `debug_workflow.py` and running + it in debug mode (see below). + + ```{important} + Before committing, confirm that `error_check_rabies.py --complete` finishes + with no errors. + ``` +5. **Commit, push and open a pull request** from your fork to the original + repository. + +## Debug interactively with Spyder + +1. Open `debug_workflow.py` in Spyder. +2. Find the scripts in your local installation and add breakpoints. Running + `import rabies; os.path.abspath(rabies.__file__)` gives you the path to the + installed package's `__init__.py`, and from there you can locate the file of + interest. +3. Run `debug_workflow.py` in debug mode until it reaches the breakpoint. + +## Add a new module to the pipeline + +RABIES workflows are structured with [Nipype](https://nipype.readthedocs.io/en/latest/) +(see also the [Nipype tutorial](https://miykael.github.io/nipype_tutorial/)). + +### Write the function as a Nipype interface + +```python +from nipype.interfaces.base import ( + traits, TraitedSpec, BaseInterfaceInputSpec, + File, BaseInterface +) + + +class NewInterfaceInputSpec(BaseInterfaceInputSpec): + # select an appropriate input type with traits.type (Dict, File, Int, ...) + input_str = traits.Str(exists=True, mandatory=True, + desc="An input string.") + + +class NewInterfaceOutputSpec(TraitedSpec): + out_file = File( + exists=True, desc="An output file.") + + +class NewInterface(BaseInterface): + """ + Describe your module. + """ + + input_spec = NewInterfaceInputSpec + output_spec = NewInterfaceOutputSpec + + def _run_interface(self, runtime): + input_str = self.inputs.input_str + + ''' + YOUR CODE + ''' + + setattr(self, 'out_file', out_file) + + return runtime + + def _list_outputs(self): + return {'out_file': getattr(self, 'out_file')} +``` + +### Wrap it in a node + +```python +from .other_script import NewInterface # if the interface is in a different script +from nipype.pipeline import engine as pe + +new_interface_node = pe.Node(NewInterface(), + name='new_interface') +``` + +A node can also be built from any plain Python function, without writing an +interface: + +```python +from nipype.pipeline import engine as pe +from nipype.interfaces.utility import Function + +new_function_node = pe.Node(Function(input_names=['input_1', 'input_2'], + output_names=['output_1', 'output_2'], + function=NewFunction), + name='new_function') +``` + +### Connect the node into a workflow + +Once the node carries out the operation you want, integrate it by linking its +inputs and outputs to other nodes. Here is a complete minimal workflow, the one +that performs slice timing correction: + +```{code-block} python +:caption: A Nipype workflow conducting slice timing correction +:linenos: + +from nipype.pipeline import engine as pe +from nipype.interfaces.utility import Function +from nipype.interfaces import utility as niu + + +def init_bold_stc_wf(name='bold_stc_wf'): + + workflow = pe.Workflow(name=name) + # intermediate node storing the workflow inputs + inputnode = pe.Node(niu.IdentityInterface( + fields=['bold_file']), name='inputnode') + # intermediate node storing the workflow outputs + outputnode = pe.Node(niu.IdentityInterface( + fields=['stc_file']), name='outputnode') + + slice_timing_correction_node = pe.Node( + Function(input_names=['in_file', 'tr', 'tpattern', 'stc_axis', + 'interp_method', 'rabies_data_type'], + output_names=['out_file'], + function=slice_timing_correction), + name='slice_timing_correction', mem_gb=1.5 * opts.scale_min_memory) + + # feed the inputnode into the STC node, and STC outputs into the outputnode + workflow.connect([ + (inputnode, slice_timing_correction_node, [('bold_file', 'in_file')]), + (slice_timing_correction_node, + outputnode, [('out_file', 'stc_file')]), + ]) + return workflow +``` + +Most contributions integrate a new interface into a pre-existing workflow +rather than creating a new one, in which case the work is connecting the right +nodes to your interface. + +```{seealso} +[Workflow reference](../reference/workflows.md) documents the existing +workflows and links to their source. +``` diff --git a/docs/how_to/index.md b/docs/how_to/index.md new file mode 100644 index 00000000..6f4fd789 --- /dev/null +++ b/docs/how_to/index.md @@ -0,0 +1,47 @@ +# How-to guides + +These guides answer "how do I ...?" for a specific goal. They assume you +already know what you want to achieve and are looking for the steps to get +there. If you are new to RABIES, work through the +[tutorial](../tutorials/index.md) first. + +## Getting RABIES running + +```{toctree} +--- +maxdepth: 1 +--- +install +run_with_containers +select_scans +change_template +``` + +## Getting good results out of it + +```{toctree} +--- +maxdepth: 1 +--- +check_orientation +troubleshoot_registration +assess_data_quality +optimise_confound_correction +minimal_preprocessing +``` + +## Contributing + +```{toctree} +--- +maxdepth: 1 +--- +contribute +``` + +```{seealso} +If your question is not covered here, ask on the +[discussion board](https://github.com/CoBrALab/RABIES/discussions), or report a +problem through the +[issue tracker](https://github.com/CoBrALab/RABIES/issues). +``` diff --git a/docs/how_to/install.md b/docs/how_to/install.md new file mode 100644 index 00000000..291eaafc --- /dev/null +++ b/docs/how_to/install.md @@ -0,0 +1,90 @@ +# How to install RABIES + +Pick the installation route that matches where you will run RABIES. + +::::{grid} 1 1 3 3 +:gutter: 2 + +:::{grid-item-card} Container +**Recommended.** Apptainer on Linux and HPC clusters, Docker elsewhere. All +dependencies included. +::: + +:::{grid-item-card} PyPI +Python package only. You install the non-Python dependencies yourself. +::: + +:::{grid-item-card} Neurodesk +Browser-based neuroimaging environment. Nothing to install locally. +::: + +:::: + +## Install as a container + +Containers package the entire computing environment, so you do not install +dependencies by hand and the software behaves identically wherever you run it. +[Apptainer](https://apptainer.org/) is generally preferred over +[Docker](https://www.docker.com) because it does not require root permissions, +which makes it usable on high performance computing clusters. + +Install Apptainer or Docker first (Apptainer publishes +[quick start guidelines](https://apptainer.org/docs/user/main/quick_start.html)), +then pull the [RABIES image](https://github.com/CoBrALab/RABIES/pkgs/container/rabies): + +::::{tab-set} + +:::{tab-item} Apptainer +:sync: apptainer + +```sh +apptainer build rabies-latest.sif docker://ghcr.io/cobralab/rabies:latest +``` + +This produces a single `.sif` file containing the whole environment. +::: + +:::{tab-item} Docker +:sync: docker + +```sh +docker pull ghcr.io/cobralab/rabies:latest +``` +::: + +:::: + +To pin a version, replace `latest` with a tag from the +[list of published images](https://github.com/CoBrALab/RABIES/pkgs/container/rabies). + +```{note} +Versions prior to 0.5.0 are not on the GitHub container registry. They remain +available on [Docker Hub](https://hub.docker.com/r/gabdesgreg/rabies). +``` + +For the execution syntax once the image is built, see +[How to handle container syntax](run_with_containers.md). + +## Install from PyPI + +RABIES is published on [PyPI](https://pypi.org/project/rabies/): + +```sh +pip install rabies +``` + +```{warning} +`pip install` gives you the Python package only. RABIES also calls out to +external neuroimaging tools, which are listed in +[`dependencies.txt`](https://github.com/CoBrALab/RABIES/blob/master/dependencies.txt) +and which you must install yourself. If you are not prepared to manage those, +use a container. +``` + +## Use RABIES on Neurodesk + +RABIES is one of the [built-in tools](https://neurodesk.github.io/applications/) +on the [Neurodesk platform](https://neurodesk.github.io/), a browser-based +neuroimaging computing environment with community-maintained prebuilt tools. +Nothing is installed on your own machine. See the +[Neurodesk documentation](https://neurodesk.github.io/docs/) to get started. diff --git a/docs/how_to/minimal_preprocessing.md b/docs/how_to/minimal_preprocessing.md new file mode 100644 index 00000000..a3e14d79 --- /dev/null +++ b/docs/how_to/minimal_preprocessing.md @@ -0,0 +1,94 @@ +# How to use RABIES on already-preprocessed data + +If your functional data was preprocessed with your own pipeline and you only +want RABIES for confound correction and analysis, you cannot skip the +preprocessing stage — but you can run it with almost everything turned off. + +Because running `rabies confound_correction --read_datasink` expects the full range of files +produced by the `preprocess` stage, it is most convenient to produce the required +intermediary outputs by running a *SHAM* preprocessing run that minimally modifies +the input data. + +## Run a SHAM preprocessing + +```sh +rabies preprocess bids_inputs/ preprocess_outputs/ \ + --anat_inho_cor method=disable \ + --bold_inho_cor method=disable \ + --commonspace_reg template_registration=no_reg,fast_commonspace=true \ + --bold2anat_coreg registration=no_reg \ + --no_HMC +``` + +```{list-table} +:header-rows: 1 +:widths: 40 60 + +* - Parameter + - What it turns off +* - `--anat_inho_cor method=disable` + - inhomogeneity correction of the structural images +* - `--bold_inho_cor method=disable` + - inhomogeneity correction of the EPI images +* - `--commonspace_reg template_registration=no_reg,fast_commonspace=true` + - unbiased template generation, and the registration to the reference atlas (an identity transform is used instead) +* - `--bold2anat_coreg registration=no_reg` + - EPI-to-structural coregistration, i.e. the susceptibility distortion correction +* - `--no_HMC` + - the *application* of head motion correction to the resampled timeseries. Head motion parameters are still estimated and remain available to `--nuisance_regressors`, `--frame_censoring` and `--data_diagnosis` +``` + +If your dataset has no structural scans, add `--bold_only`, in which case +`--anat_inho_cor` and `--bold2anat_coreg` no longer apply. + +```{warning} +`template_registration=no_reg` does not skip the resampling to commonspace — it +replaces the estimated transform with an identity transform. The commonspace +outputs are therefore only meaningful if your input images **already overlap** +with the template given to `--anat_template`, which is the file that defines +the commonspace. + +If they do not overlap, the commonspace timeseries and the atlas masks +(`--brain_mask`, `--WM_mask`, `--CSF_mask`, `--vascular_mask`) applied +downstream will not correspond to your data. +``` + +## What still happens + +A SHAM preprocessing is a minimal pass, not a strictly non-modifying one. + +Operations which alter the data but are off by default stay off, and should not +be added: `--apply_STC`, `--apply_despiking`, `--detect_dummy`, +`--log_transform`, `--anat_autobox`, `--bold_autobox` and `--oblique2card`. + +The timeseries are still resampled onto the output grid, using the identity +transforms described above. Use `--commonspace_resampling` and +`--anatomical_resampling` to control the output voxel dimensions, and +`--interpolation` to select the interpolator. + +Two further operations are applied unconditionally and cannot be turned off: + +**Negative values are clipped to zero** +: This happens when the preprocessed timeseries are written out. If your data + legitimately contains negative values — because it was already demeaned or + detrended by your own pipeline — those voxels will be set to zero. Bring in + data on a positive scale, and leave centring to `--detrending` at the + confound correction stage. + +**The output is cast** to the type given by `--data_type` +: `float32` by default. + +## Alternative: `--read_datasink` + +`rabies confound_correction --read_datasink` reads the preprocessing outputs +from the datasink folders rather than from the saved workflow graph, which +removes the need for the `.pkl` file. This requires reproducing the RABIES +[output structure and file naming](../reference/outputs.md) exactly, and is +generally more work than running a SHAM preprocessing. + +## Still not covered? + +If your use case needs settings that are not exposed, open a +[discussion](https://github.com/CoBrALab/RABIES/discussions) describing +explicitly what you need. Providing example data lets us work out an +implementation supporting your use case. diff --git a/docs/how_to/optimise_confound_correction.md b/docs/how_to/optimise_confound_correction.md new file mode 100644 index 00000000..29dace22 --- /dev/null +++ b/docs/how_to/optimise_confound_correction.md @@ -0,0 +1,64 @@ +(optim_CR)= + +# How to optimise your confound correction strategy + +There is no agreed-upon single confound correction strategy for fMRI +functional connectivity analysis. Designing a confound correction +pipeline is navigating a central trade-off: too little signal removal +risks leaving uncorrected artefacts that corrupt downstream analyses, +excessive signal removal risks removing signal of interest relating to +network activity. The ideal strategy will depend on the extent and +nature of artefacts present in a given dataset. + +Here is described a protocol introduced by {cite}`Desrosiers-Gregoire2024-ou` +to navigate these decisions on a per-dataset basis. It assumes you can +already generate and read the data quality reports — if not, start with +[How to assess data quality](assess_data_quality.md). + +## The protocol + +1. **Start with a minimal correction** and generate the data quality reports at + the analysis stage. A reasonable minimum can be frame censoring on framewise + displacement, regression of the 6 motion parameters, and spatial smoothing: + + ```sh + rabies -p MultiProc confound_correction preprocess_outputs/ confound_correction_outputs/ \ + --frame_censoring FD_censoring=true,FD_threshold=0.05 \ + --nuisance_regressors mot_6 \ + --smoothing_filter 0.3 + ``` + + ```{tip} + Starting minimal is ideal, as excessive + correction removes network activity along with the confounds, and + over-correction is harder to detect after the fact than under-correction. + ``` + +2. **Evaluate the reports**, following + [How to assess data quality](assess_data_quality.md). + +3. **Choose one additional correction**, using the table below to match what + you observed to the correction that addresses it. + +4. **Re-run confound correction with that one correction added**, regenerate + the reports, and compare. Keep the addition only if it improved the quality + outcomes. Adding one correction at a time is what makes its effect + attributable. + +5. **Repeat steps 3 and 4** until the quality outcomes are acceptable or you + have run out of options. + +## Matching observations to corrections + +```{figure} ../pics/CR_optimization_table.svg +:alt: Table relating data quality observations to the corresponding confound correction options +:width: 100% + +Guidance for prioritising additional corrections based on observations from the +data quality reports. +``` + +```{seealso} +[The confound correction workflow](confound_pipeline_target) describes every +correction step available and the order in which RABIES applies them. +``` diff --git a/docs/how_to/run_with_containers.md b/docs/how_to/run_with_containers.md new file mode 100644 index 00000000..32463a0f --- /dev/null +++ b/docs/how_to/run_with_containers.md @@ -0,0 +1,116 @@ +# How to handle container syntax + +A container has its own filesystem and cannot see your data unless you say so. +Running RABIES in a container is therefore the same as running it natively, +plus one rule: **every directory RABIES needs must be bound to a path inside +the container, and the RABIES arguments must use the container-side paths.** + +Bind directories with `-B` for Apptainer and `-v` for Docker. Both take +`host_path:container_path`, with an optional `:ro` to make the bind read-only. + +```{important} +Bind the same directories at the same container-side paths for **all three +stages**. Each stage reads the file paths recorded by the previous one, so a +path used during preprocessing must still resolve during confound correction +and analysis. Changing or dropping a bind between stages produces missing-file +errors. +``` + +## Apptainer + +### Preprocessing + +```sh +apptainer run -B $PWD/input_BIDS:/input_BIDS:ro \ + -B $PWD/preprocess_outputs:/preprocess_outputs/ \ + /path_to_apptainer_image/rabies.sif \ + -p MultiProc preprocess /input_BIDS/ /preprocess_outputs/ \ + --apply_STC --TR 1.2 \ + --commonspace_reg masking=true,brain_extraction=false,template_registration=SyN,fast_commonspace=false +``` + +`apptainer run /path_to_apptainer_image/rabies.sif` executes the image; every +argument after it is passed to RABIES and follows the ordinary +[command line syntax](../reference/cli.md). The two binds are what make the +data reachable: + +`-B $PWD/input_BIDS:/input_BIDS:ro` +: Maps your BIDS folder to `/input_BIDS` inside the container, which is why the + RABIES argument reads `/input_BIDS/`. `:ro` grants read-only access, so the + container cannot modify your raw data. + +`-B $PWD/preprocess_outputs:/preprocess_outputs/` +: Maps the desired output directory. There is no `:ro`, so the container can + write here. + +### Confound correction + +```sh +apptainer run -B $PWD/input_BIDS:/input_BIDS:ro \ + -B $PWD/preprocess_outputs:/preprocess_outputs/ \ + -B $PWD/confound_correction_outputs:/confound_correction_outputs/ \ + /path_to_apptainer_image/rabies.sif \ + -p MultiProc confound_correction /preprocess_outputs/ /confound_correction_outputs/ \ + --nuisance_regressors WM_signal CSF_signal vascular_signal mot_6 \ + --smoothing_filter 0.3 +``` + +`/input_BIDS` is still bound even though it does not appear in the RABIES +arguments — this is the rule stated above. + +### Analysis + +```sh +apptainer run -B $PWD/input_BIDS:/input_BIDS:ro \ + -B $PWD/preprocess_outputs:/preprocess_outputs/ \ + -B $PWD/confound_correction_outputs:/confound_correction_outputs/ \ + -B $PWD/analysis_outputs:/analysis_outputs/ \ + /path_to_apptainer_image/rabies.sif \ + -p MultiProc analysis /confound_correction_outputs /analysis_outputs/ \ + --group_ica apply=true,dim=30,random_seed=1 +``` + +## Docker + +The syntax mirrors Apptainer, with `-v` in place of `-B` and a few extra flags: + +```sh +docker run -it --rm --user $(id -u) \ + -v $PWD/input_BIDS:/input_BIDS:ro \ + -v $PWD/preprocess_outputs:/preprocess_outputs/ \ + ghcr.io/cobralab/rabies:latest \ + -p MultiProc preprocess /input_BIDS/ /preprocess_outputs/ \ + --apply_STC --TR 1.2 \ + --commonspace_reg masking=true,brain_extraction=false,template_registration=SyN,fast_commonspace=false +``` + +`--user $(id -u)` +: Runs as your own user id, so output files are owned by you. Without it, + Docker writes files as root and you may be unable to delete them. + +`--rm` +: Removes the container when the run finishes. + +Replace `latest` with a specific version tag for reproducible runs. + +## Using a custom atlas or seed files + +Template files, masks and seeds passed with `--anat_template`, `--brain_mask`, +`--WM_mask`, `--CSF_mask`, `--vascular_mask`, `--prior_maps` or `--seed_list` +live outside your input and output directories, so they need binds of their +own: + +```sh +apptainer run -B $PWD/input_BIDS:/input_BIDS:ro \ + -B $PWD/preprocess_outputs:/preprocess_outputs/ \ + -B $PWD/my_atlas:/atlas:ro \ + /path_to_apptainer_image/rabies.sif \ + -p MultiProc preprocess /input_BIDS/ /preprocess_outputs/ \ + --anat_template /atlas/template.nii.gz \ + --brain_mask /atlas/brain_mask.nii.gz +``` + +```{seealso} +[CoBrALab recommendations](https://github.com/CoBrALab/documentation/wiki/Running-RABIES-on-niagara) +for running RABIES on Compute Canada clusters. +``` diff --git a/docs/how_to/select_scans.md b/docs/how_to/select_scans.md new file mode 100644 index 00000000..97674448 --- /dev/null +++ b/docs/how_to/select_scans.md @@ -0,0 +1,95 @@ +# How to select which scans get processed + +By default RABIES processes every functional scan it can find in the input BIDS +directory. There are two independent ways to narrow that down: change which +files count as functional or anatomical scans, or list the individual scans to +keep or drop. + +## Change which files RABIES recognises + +RABIES identifies functional scans by the `bold` or `cbv` suffix and structural +scans by the `T1w` or `T2w` suffix. Files matching neither are ignored. + +Override this with `--bids_filter`, which takes a JSON file with the functional +filter under `func` and the structural filter under `anat`. This JSON reproduces +the defaults: + +```{code-block} json +:caption: bids_filter.json + +{ + "func": { + "suffix": ["bold", "cbv"] + }, + "anat": { + "suffix": ["T1w", "T2w"] + } +} +``` + +Pass it at the preprocessing stage: + +```sh +rabies preprocess input_BIDS/ preprocess_outputs/ --bids_filter bids_filter.json +``` + +Add any BIDS entity to either filter to be more specific. To use only the RARE +anatomical acquisition and only the resting-state functional runs: + +```{code-block} json +:caption: bids_filter.json + +{ + "func": { + "suffix": ["bold"], + "task": ["rest"] + }, + "anat": { + "suffix": ["T2w"], + "acquisition": ["RARE"] + } +} +``` + +```{tip} +Reach for `--bids_filter` when your dataset contains scans that are not meant +for this pipeline at all — a second anatomical modality, a task run alongside +the resting-state runs. Reach for `--inclusion_ids` below when you want a +subset of otherwise-eligible scans. +``` + +## Include or exclude individual scans + +`--inclusion_ids` and `--exclusion_ids` take the full paths of BOLD files. +They are execution options, so they go **before** the processing stage name, +and they can be given at any stage. + +Process only two scans: + +```sh +rabies --inclusion_ids input_BIDS/sub-001/ses-1/func/sub-001_ses-1_task-rest_bold.nii.gz \ + input_BIDS/sub-002/ses-1/func/sub-002_ses-1_task-rest_bold.nii.gz \ + -p MultiProc preprocess input_BIDS/ preprocess_outputs/ +``` + +Process everything except one scan: + +```sh +rabies --exclusion_ids input_BIDS/sub-003/ses-1/func/sub-003_ses-1_task-rest_bold.nii.gz \ + -p MultiProc preprocess input_BIDS/ preprocess_outputs/ +``` + +For longer lists, put one file path per row in a text file and pass the file: + +```sh +rabies --inclusion_ids scans_to_process.txt -p MultiProc preprocess input_BIDS/ preprocess_outputs/ +``` + +```{warning} +Do not put `--inclusion_ids` or `--exclusion_ids` immediately before the +processing stage name. The stage name gets swallowed into the list and +argument parsing fails. Put another option, such as `-p` or `--verbose`, +between them — as in every example above. + +`--inclusion_ids` and `--exclusion_ids` cannot be used together. +``` diff --git a/docs/how_to/troubleshoot_registration.md b/docs/how_to/troubleshoot_registration.md new file mode 100644 index 00000000..4d7c3d5a --- /dev/null +++ b/docs/how_to/troubleshoot_registration.md @@ -0,0 +1,89 @@ +# How to troubleshoot registration + +The RABIES defaults modify the images as little as possible and lean on the +quality of the images at acquisition. They are the right starting point, but +they do not produce a robust workflow for every dataset, and reaching good +outcomes on your own data will often mean tuning parameters. + +This guide maps failures visible in the +[preprocessing QC report](../reference/qc_outputs.md) onto the parameters that +address them. Work through it in order: inhomogeneity correction feeds masking, +and masking feeds every registration that follows, so a problem there will +reappear downstream no matter what you change at the registration step. + +```{important} +Before changing any parameter, confirm your images are correctly oriented — +see [How to check image orientation](check_orientation.md). Mis-oriented +images produce registration failures that no amount of parameter tuning will +fix. +``` + +## Inhomogeneity correction + +Relevant parameters: `--anat_inho_cor`, `--bold_inho_cor`, +`--anat_robust_inho_cor`, `--bold_robust_inho_cor`. Inspect the +`anat_inho_cor/` and `bold_inho_cor/` folders of the QC report. + +Only a subset of scans have failed masking, or the mask is partially misregistered +: Use `--anat_robust_inho_cor`/`--bold_robust_inho_cor`. These register all + corrected images together to generate a temporary template representing the + average of all scans, mask that template, and use it as the masking target in + a second pass of inhomogeneity correction, which is a more robust target. The + parameters for these options are the same as for `--commonspace_reg` below. + +Inhomogeneity biases are not completely corrected, with signal drops remaining +: Apply `multiotsu=true`. This performs a staged correction, correcting low + intensities first and iteratively including higher ones, which better handles + images with strong inhomogeneity gradients and very low intensities. + +Tissue outside the brain is causing registration failures +: If the initial correction enhanced the intensity of tissue outside the brain + and masking then fails, use `--anat_autobox`/`--bold_autobox` to crop out the + extra tissue automatically. You can also modify `otsu_thresh`, which sets the + threshold for the automatic masking during the initial correction, to select + a value more specific to brain tissue. + +A large proportion of masking failures remain — mismatched brain sizes, non-linear warps, or the mask falling outside the brain +: Apply a less stringent registration `method`, stepping down through + `SyN` → `Affine` → `Rigid` → `no_reg`. If you reach `no_reg`, you may also + have to adjust `otsu_thresh` to obtain an automatically generated brain mask + covering only brain tissue. + +## Commonspace registration and susceptibility distortion correction + +Relevant parameters: `--commonspace_reg`, `--bold2anat_coreg`. Inspect the +`commonspace_reg_wf.Anat2Unbiased/` and `commonspace_reg_wf.Unbiased2Atlas/` +folders of the QC report (or `commonspace_reg_wf.Anat2Atlas/` if you used +`fast_commonspace=true`), together with `EPI2Anat/`. + +Many scans are misregistered, or brain edges are not well matched +: First inspect the quality of inhomogeneity correction for those scans, and + follow the guidance above if the correction or brain masking was poor. If + good quality masks were obtained during inhomogeneity correction, bring them + into the registration with `masking=true`. If registration errors persist, + particularly at the brain edges, `brain_extraction=true` and `keep_mask_after_extract=true` further constrain + edge matching by removing tissue outside the brain. + + ```{warning} + The quality of brain edge delineation depends on the masks derived during + inhomogeneity correction, so `brain_extraction=true` is only as good as that + earlier step. Fix masking first. + ``` + +Scans have incomplete brain coverage, and surrounding tissue is stretched to fill the gap +: Non-linear registration assumes corresponding anatomy between the moving + image and the target. When brain regions are missing — the cerebellum or + olfactory bulbs are the usual cases — the surrounding tissue may be + improperly stretched to fill the missing area. `brain_extraction=true,keep_mask_after_extract=true` can + largely mitigate this. + +## When registration cannot be salvaged + +If a scan cannot be registered acceptably, exclude it rather than analysing it: +see [How to select which scans get processed](select_scans.md). + +```{seealso} +- [Preprocessing QC outputs](../reference/qc_outputs.md) — what each QC folder shows +- [The preprocessing workflow](../explanation/preprocessing.md) — what each registration step is for +- [`preprocess` options](../reference/cli.md) — the complete parameter list, with all accepted values +``` diff --git a/docs/index.md b/docs/index.md index 494d2cda..abfb2979 100644 --- a/docs/index.md +++ b/docs/index.md @@ -1,21 +1,81 @@ ```{include} ../README.md ``` +## Finding your way around this documentation + +This documentation is organised around what you are trying to do right now. + +::::{grid} 1 1 2 2 +:gutter: 3 + +:::{grid-item-card} {octicon}`rocket;1.5em;sd-mr-1` Tutorial +:link: tutorials/index +:link-type: doc + +**Start here if RABIES is new to you.** A guided run of the complete +pipeline on a small example dataset, from raw BIDS input to a +connectivity map. +::: + +:::{grid-item-card} {octicon}`tools;1.5em;sd-mr-1` How-to guides +:link: how_to/index +:link-type: doc + +**Practical recipes for a specific goal.** Installing RABIES, running it +in a container, tuning a failed registration, designing a confound +correction strategy, contributing code. +::: + +:::{grid-item-card} {octicon}`book;1.5em;sd-mr-1` Reference +:link: reference/index +:link-type: doc + +**Look something up.** Command line options, input requirements, every +output file, and the precise definition of every metric. +::: + +:::{grid-item-card} {octicon}`light-bulb;1.5em;sd-mr-1` Explanation +:link: explanation/index +:link-type: doc + +**Understand how RABIES works and why.** The preprocessing and confound +correction workflows, the analyses, and the data quality framework. +::: + +:::: + ```{toctree} --- -maxdepth: 3 -caption: Content ---- -installation.md -running_the_software.md -preprocessing.md -preproc_QC.md -confound_correction.md -analysis.md -analysis_QC.md -outputs.md -metrics.md -faq.md -troubleshooting.md -contributing.md -bibliography.md +maxdepth: 2 +hidden: +caption: Getting started +--- +tutorials/index +``` + +```{toctree} +--- +maxdepth: 2 +hidden: +caption: How-to guides +--- +how_to/index +``` + +```{toctree} +--- +maxdepth: 2 +hidden: +caption: Reference +--- +reference/index +``` + +```{toctree} +--- +maxdepth: 2 +hidden: +caption: Explanation +--- +explanation/index +``` diff --git a/docs/installation.md b/docs/installation.md deleted file mode 100644 index 3eb3ea94..00000000 --- a/docs/installation.md +++ /dev/null @@ -1,25 +0,0 @@ -# Installation - -## Container (Apptainer/Docker) \*\*RECOMMENDED\*\* -For most uses, we recommend instead using a containerized installation with [Apptainer](https://apptainer.org/) when possible on a Linux system (here's their [quick start guidelines](https://apptainer.org/docs/user/main/quick_start.html)), or [Docker](https://www.docker.com) on other platforms. Containers allow to build entire computing environments, grouping all dependencies required to run the software. This in turn reduces the burden of installing dependencies manually and ensures reproducible behavior of the software. Apptainer is generally preferred over Docker since root permissions are not required, and is thus generally compatible across computing platforms (e.g. high performance computing clusters). - -A [containerized version](https://github.com/CoBrALab/RABIES/pkgs/container/rabies) of RABIES is available from Github. After installing Apptainer or Docker, the following command will pull and build the container: -* Install Apptainer .sif file: -``` -apptainer build rabies-latest.sif docker://ghcr.io/cobralab/rabies:latest -``` -* Install Docker image: -``` -docker pull ghcr.io/cobralab/rabies:latest -``` -A specific tag version can be selected (instead of `latest`) from the [list online](https://github.com/CoBrALab/RABIES/pkgs/container/rabies). Versions prior to 0.5.0 are found on [Docker Hub](https://hub.docker.com/r/gabdesgreg/rabies). - -## PyPi -The software is available on [PyPi](https://pypi.org/project/rabies/), which makes the rabies python package widely accessible with -``` -pip install rabies -``` -However, this does not account for non-python dependencies found in `dependencies.txt`. - -## Neurodesk -RABIES is also made available on the [Neurodesk platform](https://neurodesk.github.io/), as part of the [built-in tools](https://neurodesk.github.io/applications/) for neuroimaging. The Neurodesk platform allows for an entirely browser-based neuroimaging computing environment, with pre-built neuroimaging tools from the community, and aims at reducing needs for manual development of computing environments and at improving reproducible neuroimaging. More details on Neurodesk here . \ No newline at end of file diff --git a/docs/metrics.md b/docs/metrics.md deleted file mode 100644 index 4fe417a8..00000000 --- a/docs/metrics.md +++ /dev/null @@ -1,71 +0,0 @@ -# Metric definitions -(metrics_target)= -On this page, the root-mean square (RMS) corresponds to $||x||_2 = \sqrt{\frac{1}{n}\sum_{i=1}^{n}x_i^2}$ - -(regressor_target)= -## Nuisance regressors for confound regression -* **mot_6**: Corresponds to the head motion translation and rotation parameters. Prior to the regression, the motion regressors are also subjected to the same frame censoring, detrending and frequency filtering which were applied to the BOLD timeseries to avoid the re-introduction of previously corrected confounds, as recommend in {cite}`Power2014-yf` and {cite}`Lindquist2019-lq`. -* **mot_24**: Corresponds to the 6 motion parameters together with their temporal derivatives, and 12 additional parameters are obtained by taking the squared terms (i.e. Friston 24 parameters {cite}`Friston1996-sa`) -$$ -mot24_t = [mot6_t,(mot6_t-mot6_{t-1}),(mot6_t)^2,(mot6_t-mot6_{t-1})^2] -$$ -with $mot24_t$ representing the list of 24 regressors for timepoint $t$. As with mot_6, the 24 regressors are additionally subjected to censoring, detrending and frequency filtering if applied on BOLD. -* **WM/CSF/vascular/global signal**: The mean signal is computed within the corresponding brain mask (WM,CSF,vascular or whole-brain mask) from the partially cleaned timeseries (i.e. after the confound correction steps 1-4 up to frequency filtering). -* **aCompCor_percent**: Principal component timecourses are derived from timeseries within the combined WM and CSF masks (aCompCor technique {cite}`Muschelli2014-vi`). From the timeseries within the WM/CSF masks $Y_{WM/CSF}$, a principal component analysis (PCA) decomposition is conducted to derive -$$ -Y_{WM/CSF} = W_{aCompCor}C^T -$$ -with $C$ corresponding to a set of spatial principal components, and $W$ to their associated loadings across time. The set of first components explaining 50% of the variance are kept, and their loadings $W_{aCompCor}$ provide the set of aCompCor nuisance regressors. PCA is conducted on the partially cleaned timeseries (i.e. after the confound correction steps 1-4 up to frequency filtering). -* **aCompCor_5**: Same as **aCompCor_percent**, but the first 5 components are kept instead of a set explaining 50% of the variance. - - -## Temporal scan diagnosis - -(mot6_target)= -* **Head motion translation and rotation parameters**: Corresponds to 3 rotations (Euler angles in radians) and 3 translations (in mm) measured for head motion realignment at each timeframe. -(FD_target)= -* **Framewise displacement**: For each timepoint, corresponds to the displacement (mean across the brain voxels) between the current and the next frame. For each brain voxel within the referential space for head realignment (i.e. the [3D EPI](3D_EPI_target) which was provided as reference for realignment) and for each timepoint, the inverse transform of the head motion parameters (from the corresponding timepoint) is applied to obtain the voxel position pre-motion correction. Framewise displacement can then be computed for each voxel by computing the Euclidean distance between the positions pre-motion correction for the current and next timepoints. Thus, the mean framewise displacement $FD_t$ at timepoint $t$ is computed as -$$ -FD_t = \frac{1}{n}\sum_{i=1}^{n}\sqrt{(x_{i,t+1}-x_{i,t})^2+(y_{i,t+1}-y_{i,t})^2+(z_{i,t+1}-z_{i,t})^2} -$$ -using the 3D $x$,$y$ and $z$ spatial coordinates in mm for timepoints $t$ and $t+1$ and for each voxel indices $i$. Framewise displacement for the last frame (which has no future timepoint) is set to 0. -(DVARS_target)= -* **DVARS**: represents the estimation of temporal shifts in global signal at each timepoint, measured as the root-mean-square of the timeseries’ temporal derivative -$$ -DVARS_t = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(Y_{i,t}-Y_{i,t-1})^2} -$$ -where $Y_{i,t}$ corresponds to the BOLD signal in brain voxel $i$ at timepoint $t$. The first timepoint is set to 0 (has no previous timepoint). -* **Framewise distance from mean**: This corresponds to the mean square error (MSE) between each frame and the average EPI, with the average computed as the tri-mean across time voxelwise. The input to this function is the resampled timeseries outputed from the `preprocess` stage. -* **Whole-brain/Edge/WM/CSF mask**: The mean signal across a certain brain mask, either whole-brain (i.e. the global signal), WM, CSF or brain edge. -* **$CR_{var}$**: The variance estimated by confound regression is computed for each timepoint. This is done by taking $CR_{var} = RMS(Y_{CR})$ across voxels at each timepoints, where $Y_{CR}$ is the [predicted confound timeseries](CR_target). -* **CR $R^2$**: Represents the proportion of variance explained (and removed) by confound regression. This is obtained with $CR_{R^2}= 1-\frac{Var(\hat{Y})}{Var(Y)}$ at each timepoint, where $Y$ and $\hat{Y}$ are the timeseries pre- and post-regression respectively, and $Var(x) = \frac{1}{n}\sum_{i=1}^{n}(x_i - \mu_x)^2$ calculates the variance, with $\mu$ as the mean. -* **Mean amplitude**: A set of timecourse are averaged as $\frac{1}{n}\sum_{i=1}^{n}|X_i|$, where $X_i$ is the timecourse $i$. Timecourses can correspond to either of the following sets: - * DR confounds: timecourses from the first stage of dual regression, using confound components provided to `--prior_confound_idx`. - * DR networks: network timecourses from the first stage of dual regression as specified with `--prior_bold_idx`. - * SBC networks: network timecourses derived from the set of seeds provided in `--seed_list`. - -## Spatial scan diagnosis - -* **BOLDSD**: The temporal standard deviation is computed for each voxel from the BOLD timeseries. -* **CRSD**: The temporal standard deviation computed on each voxel from the predicted confound timeseries during confound regression (i.e. [$Y_{CR}$](CR_target)). -* **CR R2**: The proportion of variance explained by confound regression at each voxel. This is obtained with $CR_{R^2}= 1-\frac{Var(\hat{Y})}{Var(Y)}$ at each voxel, where $Y$ and $\hat{Y}$ are the timeseries pre- and post-regression respectively, and $Var(x) = \frac{1}{n}\sum_{i=1}^{n}(x_i - \mu_x)^2$ is variance of $x$, with $\mu$ as the mean. -* **Global signal covariance (GScov)**: The covariance between the global signal and the timeseries at each voxel is measured as $GS_{cov} = \frac{1}{n}\sum_{t=1}^{n}Y_t \times GS_t$, where $GS_t = \frac{1}{n}\sum_{i=1}^{n}Y_i$, i.e. the mean across all brain voxels for a given timepoint. -* **DR network X**: The linear coefficients resulting from the [second regression with dual regression](DR_target), corresponding to a network amplitude map (for the Xth network specified for analysis with `--prior_bold_idx`). -* **SBC network X**: The voxelwise correlation coefficients (pearson's r) estimated with seed-based connectivity (for the Xth seed provided for analysis with `--seed_list`). - - -## Distribution plot -(dist_plot_metrics)= - -* **Network amplitude**: The overall network amplitude is summarized by computing the L2-norm across a network connectivity map from a subject-level analysis (such a map can be derived from seed-based correlation, or correspond to the linear coefficients from the [second regression ${\beta}_{SM}$](DR_target) for dual regression) -* **Network specificity**: The network map (seed-based or dual regression) and the corresponding canonical network map are thresholded to include the top X% of voxels with highest connectivity (X% is defined by `--brainmap_percent_threshold`), and the overlap of the thresholded area is computed using Dice overlap. For dual regression, the 'canonical network' map will consist of the original ICA component corresponding to that network provided with `--prior_maps`, and for seed-based connectivity, the reference network maps are provided using the `--seed_prior_list` parameter. -* **Dual regression confound correlation**: The timecourse for a single network (from a seed or dual regression) is correlated with the timecourse from each confound component (provided using `--prior_confound_idx`) modelled through dual regression, then the absolute mean correlation is computed to obtain the average amplitude of confound correlations for this specific network analysis. -* **FD-DVARS corr.**: For each scan, this QC metric corresponds to the correlation between the framewise displacement (FD) timecourse and the DVARS timecourse **computed post-confound correction** (i.e. this is not the DVARS plotted in the temporal diagnosis figure). Crucially, censored timeframes are excluded from both timecourses, and DVARS was recomputed after applying confound correction so that this metric represents *uncorrected/residual* associations between spontaneous motion and the cleaned global signal fluctuations. -* **Total $CR_{SD}$**: The total standard deviation across the [predicted confound timeseries $Y_{CR}$](CR_target). -* **Mean framewise displacement**: The mean framewise displacement computed across time (only including frames after censoring applied for confound correction). -* **Temporal degrees of freedom**: The remaining degrees of freedom post-confound correction are calculated as `tDOF = Original number of timepoints - Number of censored timepoints - Number of AROMA components removed - Number of nuisance regressors`. - -## Group statistical QC report -(group_QC_metrics)= -* **Specificity of network variability**: similarly to network specificity in the distribution plot, the network variability map and the corresponding canonical network map are thresholded to include the top X% of voxels (X% is defined by `--brainmap_percent_threshold`), and then the overlap is estimated using Dice overlap. -* **Mean confound correlation**: for each confound correlation map (either $CR_{SD}$, mean FD or tDOF), the mean is computed across voxels included within the thresholded area of the canonical network map, to obtain a mean correlation within the network's core region. \ No newline at end of file diff --git a/docs/nested_docs/distribution_plot.md b/docs/nested_docs/distribution_plot.md deleted file mode 100644 index 137aba2c..00000000 --- a/docs/nested_docs/distribution_plot.md +++ /dev/null @@ -1,16 +0,0 @@ -# Distribution plot - -(dist_plot_target)= - -![](../pics/distribution_plot.png) - -The distribution plot allows visualizing the distribution of data quality measures across the dataset, where measures of network connectivity (specificity and amplitude) are contrasted with measures of confounds across samples (each point in the plot is a scan). Data points labelled in gray were removed using `--scan_QC_thresholds`, where the gray dotted lines correspond to the QC thresholds selected for network specificity (Dice overlap) and DR confound correlation. Among the remaining samples and for each metric separately, scan presenting outlier values were detected based on a modified Z-score threshold (set with `--outlier_threshold`, 3.5 by default) and labelled in orange. The derivation of the quality metrics is described in details on the [metrics documentation](dist_plot_metrics). -
- -The report was designed to subserve two main functions: 1. Inspect that network specificity is sufficient and the temporal correlation with confounds (i.e. DR confound corr.) minimal, and set thresholds for scan inclusion using `--scan_QC_thresholds` (top right subplot, more details on this below), and 2. complement the group statistical report to visualize the association between connectivity and the three confound measures included in the report ($CR_{SD}$, mean FD and tDOF). In the later case, it can be possible for instance to determine whether a group-wise correlation in statistical report is driven by outliers. - -## Scan-level thresholds based on network specificity and confound temporal correlation - -![](../pics/scan_QC_thresholds.png) - -The measures of network specificity (using Dice overlap) and temporal correlation with confounds (where confound timecourses are extracted using confound components specified with `--conf_prior_idx` and measured through dual regression) were defined in {cite}`Desrosiers-Gregoire2024-ou` for conducting scan-level QC (the figure above is reproduced from the study). They were selected as ideal measures for quantifying issues of network detectability and spurious connectivity (the figure above demonstrate how [categories of scan quality outcomes](quality_marker_target) can be distinguished with these metrics), and applying inclusion thresholds to select scans which respect assumptions for network detectability and minimal effects from confounds. diff --git a/docs/nested_docs/group_stats.md b/docs/nested_docs/group_stats.md deleted file mode 100644 index c6f69658..00000000 --- a/docs/nested_docs/group_stats.md +++ /dev/null @@ -1,15 +0,0 @@ -# Group stats - -(group_stats_target)= - -![](../pics/group_stats_QC.svg) - -Inspecting scan-level features is insufficient to conclude that inter-scan variability in connectivity isn't itself impacted (which is of primary interest for group analysis). This final report is aimed at inspecting features of connectivity variability at the group level, and focuses on two aspects: - -1. **Specificity of network variability:** the standard deviation in connectivity across scan is computed voxelwise. This allows to visualize the spatial contrast of network variability. If primarily driven by network connectivity, the contrast should reflect the anatomical extent of the network of interest (as in the example above for the mouse somatomotor network), or otherwise may display spurious or absent features. For more details on the development of this metric, consult {cite}`Desrosiers-Gregoire2024-ou`. - - **Relationship to sample size**: {cite}`Desrosiers-Gregoire2024-ou` demonstrate that the contrast of the network variability map depends on sample size. If network connectivity is observed in individual scans, but not in this statistical report, increasing sample size may improve this contrast. -2. **Correlation with confounds:** Connectivity is correlated across subject, for each voxel, with either of the three measures of confound included: variance explained from confound correction at a given voxel ($CR_{SD}$, see [predicted confound timeseries $Y_{CR}$](CR_target)), mean framewise displacement (FD), or temporal degrees of freedom. This allows establishing the importance of the association with potential confounds. What constitute a 'concerning' correlation may depend on the study, and the effect size of interest (i.e. is the effect size of interest much higher or similar to the effect size of confounds?). - -**Quantitative CSV report**: A CSV file is also automatically generated along the figure, which records a quantitative assessment of these two aspects. More specifically, the overlap between the network variability map and the reference network map is measuring using Dice overlap, and for confound measures, the mean correlation is measured within the area of the network (consult the [metric details elsewhere](group_QC_metrics)). These measures can be referred to for a quantitative summary instead (although visualization is preferred, as the Dice overlap for network variability may not perfectly distinguish network and spurious features). - -**IMPORTANT**: the validity of this report is dependent on whether [scan-level assumptions](dist_plot_target) of network detectability and minimal confound effects are met. This is because either the lack of network activity or spurious effects in a subset of scan can drive 'apparent' network variability, since there will be differences in the presence VS absence of the network across scans, but these differences may be actually driven by data quality divergences. \ No newline at end of file diff --git a/docs/nested_docs/optim_CR.md b/docs/nested_docs/optim_CR.md deleted file mode 100644 index fd8ed221..00000000 --- a/docs/nested_docs/optim_CR.md +++ /dev/null @@ -1,15 +0,0 @@ -# Optimization of confound correction strategy - -(optim_CR)= - -On this page is a procedure for improving confound correction design based on observations from the data quality assessment reports. These recommendations were originally developed in {cite}`Desrosiers-Gregoire2024-ou`, and consist of a stepwise protocol where confound correction is improved incrementally while referring to data quality reports and the table found on this page, relating data quality features to corresponding corrections. The protocol is as follows: - -1. Initiate a **minimal** confound correction, and generate data quality reports at the analysis stage. Correction should be minimal at first to mitigate potential issues of over-correction, where network activity itself can be removed by excessive correction. A minimal correction can consist of applying frame censoring using framewise displacement and the regression of 6 motion parameters together with spatial smoothing. -2. Evaluation of the data quality reports (as described in the [guidelines on the main page](analysis_QC_target)) -3. The most sensible additional correction is selected based on the observations and using the table below. -4. The confound correction pipeline stage is re-run with **one** additional correction at a time, and the data quality reports are re-evaluated. Only a single correction is tested at a time so its impact can be evaluated, and the correction is only kept if there were beneficial impacts. -5. Repeat 3. and 4. until desirable quality outcomes are met, or no further options are left for confound correction. - -The table below offers guidance for prioritizing additional corrections based on observations from the data quality reports. The confound correction workflow and the various strategies available are described elsewhere in the [confound correction pipeline](confound_pipeline_target). - -![](../pics/CR_optimization_table.svg) \ No newline at end of file diff --git a/docs/nested_docs/registration_troubleshoot.md b/docs/nested_docs/registration_troubleshoot.md deleted file mode 100644 index 49634f32..00000000 --- a/docs/nested_docs/registration_troubleshoot.md +++ /dev/null @@ -1,17 +0,0 @@ -# Recommendations for registration troubleshooting -When first attemting preprocessing with RABIES, we recommend following the default parameters as they involve less stringent modifications of the images and mostly rely on the original quality of the MR images at acquisition. However, the default parameters do not offer a generalizable robust workflow for every datasets, and to reach ideal outcomes, the workflow parameters may require tuning. We provide below recommendations for common types of registration failures that may be found from the QC report. - - -## Inhomogeneity correction (anat or BOLD) `--anat_inho_cor`, `--bold_inho_cor`, `--anat_robust_inho_cor`, `--bold_robust_inho_cor` - -* **Only a subset of the scans have failed masking, or the mask is partially misregistered:** Consider using the `--anat_robust_inho_cor/--bold_robust_inho_cor` option, which will register all corrected images to generate a temporary template representing the average of all scans, and this template is then itself masked, and becomes the new target for masking during a second iteration of inhomogeneity correction. This should provide a more robust registration target for masking. The parameters for handling this setp are the same as `--commonspace_reg` below. -* **The inhomogeneity biases are not completely corrected:** if you observe that drops in signal are still present after the connection, you should consider applying `multiotsu=true`. This option will better correct low intensities in an image with important signal drops. -* **Tissue outside the brain is provoking registration failures:** if the intensity of tissue outside the brain was enhanced during the initial inhomogeneity correction and leads to masking failures, you can consider using `--anat_autobox/--bold_autobox` which can automatically crop out extra tissue. You can also modify the `otsu_thresh` to set the threshold for the automatic masking during the initial correction, and attempt to select a threshold that is more specific to the brain tissue. -* **There are still a large proportion of masking failures (mismatched brain sizes or non-linear wraps, or mask outside of the brain):** Consider applying a less stringent registration `method`, going down from `SyN` -> `Affine` -> `Rigid` -> `no_reg` . If `no_reg` is selected, you may have to also adjust the `otsu_thresh` to obtain an automatically-generated brain mask covering only the brain tissues. - - -## Commonspace registration `--commonspace_reg` or susceptibility distortion correction `--bold2anat_coreg` - -* **Many scans are misregistered, or brain edges are not well-matched:** First, inspect the quality of inhomogeneity correction for those scans, and refer to instructions above if the correction or brain masking was poor. If good quality masks were obtained during inhomogeneity correction, they can be used to improve registration quality by using `masking=true`. If registration errors persist, in particular if brain edges are not well-matched, `brain_extraction=true` can be used to further constrain the matching of brain edges after removing tissue outside the brain. However, the quality of brain edge delineation depends on masks derived during inhomogeneity correction, so this option depends on high quality masking during this previous step. -* **Scans have incomplete brain coverage (e.g. cerebellum/olfactory bulbs), and surrounding brain tissue is streched to fill in missing regions:** The non-linear registration assumes corresponding brain anatomy between the moving image and the target. If brain regions are missing, the surrounding tissue may be improperly stretched to fill missing areas. Using the `brain_extraction=true` can largely mitigate this issue. - diff --git a/docs/nested_docs/scan_diagnosis.md b/docs/nested_docs/scan_diagnosis.md deleted file mode 100644 index 738bb11a..00000000 --- a/docs/nested_docs/scan_diagnosis.md +++ /dev/null @@ -1,55 +0,0 @@ -# Scan diagnosis report -(diagnosis_target)= - -By executing `--data_diagnosis` at the analysis stage of the pipeline, a set of visual reports are generated to support data quality assessment in relationship to connectivity analysis. Here, the *spatiotemporal diagnosis* report is described. The diagnosis is a visual report generated for each scan independently after conducting dual regression or seed-based connectivity analysis. It will display a large set of temporal and spatial features for the scan supporting the assessment of potential data quality issues, and whether network connectivity is impacted. Unless specified otherwise, all metrics are computed from fMRI timeseries after the confound correction stage. This page first covers an example of the report with the description for the set of features, and second provides guidance for interpreting the report. - - -## Spatiotemporal diagnosis -![](../pics/spatiotemporal_diagnosis.png) -Above is an example of the report (files generated into the `data_diagnosis_datasink/figure_temporal_diagnosis/` and `data_diagnosis_datasink/figure_spatial_diagnosis/` folders) for a scan with little confound signatures and clear network connectivity. Each spatial map is represented along 6 cortical slices, overlapped onto the anatomical template in common space. The network maps from dual regression (DR) or seed-based connectivity (SBC) are thresholded to include the top X% of the voxels with the highest values (X% is defined by `--brainmap_percent_threshold`). In this example, both dual regression and seed-based connectivity was conducted, where DR network 0 and SBC network 1 correspond to analysis of the somatomotor network, whereas DR network 1 and SBC network 0 correspond to the default mode network. Below we detail the interpretation of each feature included in the diagnosis (whereas the detailed computations for each metric are further described in the [Metric definitions](metrics_target) page): - -### Temporal diagnosis from top to bottom: -* **Power spectrum:** the frequency power spectrum is displayed to assess the dominant frequency profile. This always shown averaged across the whole brain, but it is also possible to display the spectrum for seeds using `--plot_seed_frequencies`. - -* **Carpet plot:** the entire fMRI timeseries are displayed in a time by voxel 2D matrix. This allows to visualize global fluctuations in signal intensity, which can be a proxy for various global artefacts {cite}`Power2017-wn`. Censored frames are labelled in red. - -* **The translation and rotation head motion parameters:** those are the 6 rigid body parameters estimated during preprocessing, and allow tracking of head position across scan duration. - -* **Framewise distance from the mean image**: This corresponds to the mean square error (MSE) between each frame and the average EPI (computed as the tri-mean across time voxelwise), computed prior to confound correction to map changes in brain constrast over time. This is a complementary measure of motion, since that metric represents changes in raw image contrast across time that remain post-head motion realignment, and thus might not be well captured by framewise displacement. This metric is also outputted using `--hmc_qc_report` during `preprocess` stage, and it is possible to censore using this metric with `--frame_censoring` at `confound_correction` stage. - -* **The framewise displacement and the temporal shifts in global signal from the root-mean-square of the timeseries’ temporal derivative (DVARS) {cite}`Power2012-ji`:** Framewise displacement quantifies movement between consecutive frames, which reveals the timing and amplitude of spontaneous motion, whereas DVARS reveals shifts in global fMRI signal intensities (which can also indicate suspicious spikes in signal). - -* **Variance captured and removed by confound regression:** The linear model fitted during nuisance regression provides an estimate of 'nuisance timeseries', from which it is possible to quantify the variance removed at each frame. The temporal diagnosis displays the absolute variance by computing the RMS of nuisance timeseries across voxels for each frame, as well as the relative variance removed by computing the R2 at each frame. These features allow both to visualize confound effects, and evaluate whether confound regression appropriately modelled confounds detected from other temporal features. - -* **Mean signals within anatomical masks:** The mean timecourse derived within a set of anatomical masks (whole-brain, edge, white matter and CSF masks) is displayed. The edge, WM and CSF regions are susceptible to motion and physiological confounds {cite}`Pruim2015-nm`, and can be contrasted with the whole-brain global signal to assess potential artefactual effects across the brain. - -* **Mean amplitude of network VS confound timecourses:** The averaged timecourse between network analyses and confound sources are compared to assess whether network amplitude is spurious (i.e. correlated with confound timecourse). To model confound timecourses, dual regression analysis is conducted with a complete set of components from Independent Component Analysis representing a mixture of networks and confounds from various origins, and the timecourses from confound components are compiled to summarize a broad set of potential confounds (by default, RABIES [this set](https://zenodo.org/record/5118030/files/melodic_IC.nii.gz) of ICA components for mice). - -### Spatial diagnosis from top to bottom: -* **Spatial distribution in signal variability (BOLDSD):** The first spatial feature of the diagnosis is the signal variability (standard deviation) at each voxel. This map offers an index of whether significant confounds are contributing to the signal (see other examples in **Interpretation of the report and main features to inspect**). Without the influence from confounds, as in this example, signal variability is largely homogeneous. - -* **Confound regression variance explained (CRSD):** The variance explained from confound regression is quantified at each voxel by taking the standard deviation from the modelled confound timecourse. This allows to contrast spatially the amplitude of confound effects. This feature can specifically delineate the presence of confounds and identify the type of confound. In this example, minor motion signatures are identified. - -* **Confound regression variance explained proportion:** Similar to CRSD, but showing instead the proportion of variance explained (R2). - -* **Global signal covariance:** This map displays the covariance of each voxel with the global signal. The contrast from this map allows to evaluate the predominant source of global signal fluctuation, which can take various forms depending on the contributions from neural network and confounds (see examples below in **Interpretation of the report and main features to inspect**). In the ideal case, there is predominant contrast found in gray matter, with a shape reminescent of brain network, as in the example shown above. - -* **Network spatial maps:** Finally, the diagnosis shows the spatial network maps fitted using dual regression (or seed-based analysis) from the selected set of brain networks of interest (in this case the somatomotor and default mode networks). These fits provide insights into the quality of network analysis, and how they may affect downstream statistical analyses. - -* Note that $CR_{SD}$ and $CR_{R^2}$ are computed from the regressors specified with `--nuisance_regressors` during confound correction. If no regressors are specified, $CR_{SD}$ and $CR_{R^2}$ are still estimated using regression on the 6 motion parameters, but that regression is not applied to remove signal from the timeseries. - -## Interpretation of the report and main features to inspect -(quality_marker_target)= -![](../pics/diagnosis_key_markers.png) - -A subset of the features in the spatiotemporal diagnosis are most crucial in determining scan quality in relationship to connectivity analysis, and are displayed above across 4 main categories of scan quality. Below we describe the key role of these four features in relationship to those 4 scan categories: - -* **BOLD variability:** The resulting BOLD variability map presents an homogeneous contrast in uncorrupted scans, and can otherwise reveal the anatomical signature of confounds, thus allowing to identify the type of confound. - -* **Global signal covariance:** The global signal covariance map is sensitive to both non-neural confounds (e.g. the spurious category) and network signatures (e.g. the specific category). The global signal covariance thus reflects whether network or confound sources dominate coordinated fluctuations, and can delineate the most likely contributors to downstream connectivity measures. - -* **Network map:** Allows inspecting whether the expected anatomical features of the network are effectively captured (i.e. network specificity). This is most crucial in ensuring that the network is not absent (see the absent category), or to ensure that the network shape is not distorded with spurious features (see spurious category). - -* **Network and confound timecourses:** Finally, the respective timecourses for networks and confounds can be compared to reveal direct relationships between network amplitude and confounds in the temporal domain. Although this metric does not describe the type of confound, it is the most direct indicator of spurious connectivity. It is an important complement to the inspection of network shape, since spurious effects may only affect amplitude with minimal impact on shape. - -These 4 features are sufficient to capture the essential characteristics of network detectability and spurious connectivity at the single scan level. The remaining features from the spatiotemporal diagnosis provide additional details regarding timeseries properties, the motion parameters, or confound regression, and can further support characterizing the specific origin of confounds (e.g. determining that a correlation between network and confound timecourse is originating from framewise displacement (i.e. motion)). diff --git a/docs/outputs.md b/docs/outputs.md deleted file mode 100644 index 21e4597b..00000000 --- a/docs/outputs.md +++ /dev/null @@ -1,99 +0,0 @@ -# Understanding the Outputs - -In this section, there is a description for all the output files provided at each processing stage. Important outputs from RABIES are stored into `datasink/` folders, which will be generated in the output folder specified at execution. - -## Preprocessing Outputs - -Multiple datasink folders are generated during preprocessing for different output types: `anat_datasink/`, `bold_datasink/`, `unbiased_template_datasink/`, `transforms_datasink/` and `confounds_datasink/`. - -- `anat_datasink/`: Includes the inhomogeneity-correction anatomical scans. - - `anat_preproc/`: anatomical scans after inhomogeneity correction - -- `bold_datasink/`: Includes all outputs related to the functional scans, where files are either resampled onto the native or commonspace of the EPI. The native space outputs are resampled over the anatomical scan from each corresponding MRI session, whereas the commonspace outputs are resampled over the reference atlas (the original EPI voxel resolution is unchanged during resampling unless specified otherwise in the RABIES command). - - `native_bold/`: preprocessed EPI timeseries resampled to nativespace - - `native_brain_mask/`: brain mask in nativespace - - `native_WM_mask/`: WM mask in nativespace - - `native_CSF_mask/`: CSF mask in nativespace - - `native_labels/`: atlas labels in nativespace - - `native_bold_ref/`: a volumetric 3D EPI average generated from the 4D `native_bold/` - - `commonspace_bold/`: preprocessed EPI timeseries resampled to commonspace - - `commonspace_mask/`: brain mask in commonspace - - `commonspace_WM_mask/`: WM mask in commonspace - - `commonspace_CSF_mask/`: CSF mask in commonspace - - `commonspace_vascular_mask/`: vascular mask in commonspace - - `commonspace_labels/`: atlas labels in commonspace - - `commonspace_resampled_template/`: the commonspace anatomical template, resampled to the EPI's dimensions - - `input_bold/`: the raw EPI scans provided as inputs in the BIDS data folder - - `initial_bold_ref/`: the initial volumetric 3D EPI average generated from the 4D `input_bold/` - - `raw_brain_mask/`: brain mask resampled onto the 4D `input_bold/` - - `inho_cor_bold/`: the volumetric 3D EPI (`initial_bold_ref/`) after inhomogeneity correction, which is later used for registration of the EPI - - `inho_cor_bold_warped2anat/`: inho_cor_bold after co-registration to the associated anatomical image (`anat_preproc/`) - - `std_map_preprocess/`: the temporal standard deviation at each voxel on the `commonspace_bold/` - - `tSNR_map_preprocess/`: the temporal signal-to-noise ratio (tSNR) of the `commonspace_bold/` - -- `unbiased_template_datasink/`: Outputs related to the generation of the unbiased template using https://github.com/CoBrALab/optimized_antsMultivariateTemplateConstruction. The unbiased template corresponds to the average of all anatomical (or functional with `--bold_only`) scans after their alignment. - - `unbiased_template/`: the unbiased template generated from the input dataset scans - - `warped_unbiased_template/`: the unbiased template, registered to the reference atlas in commonspace - -- `transforms_datasink/`: datasink for all the relevant transform files resampling between the different spaces. The bold_to_anat registration transformed the raw EPI to overlap with the anatomical image, correcting for susceptibility distortions, which corresponds to the native space. The native_to_unbiased registration overlaps every scans to the generated unbiased template, and then the unbiased_to_atlas corresponds to the registration of the unbiased template with the reference atlas, which defines the commonspace. - - `bold_to_anat_affine/`: affine transforms from the EPI co-registration to the anatomical image - - `bold_to_anat_warp/`: non-linear transforms from the EPI co-registration to the anatomical image - - `bold_to_anat_inverse_warp/`: inverse of the non-linear transforms from the EPI co-registration to the anatomical image - - `native_to_unbiased_affine/`: affine transforms for the alignment between native space and the unbiased template - - `native_to_unbiased_warp/`: non-linear transforms for the alignment between native space and the unbiased template - - `native_to_unbiased_inverse_warp/`: inverse of the non-linear transforms for the alignment between native space and the unbiased template - - `unbiased_to_atlas_affine/`: affine transforms for the alignment between unbiased template and the atlas in commonspace - - `unbiased_to_atlas_warp/`: non-linear transforms for the alignment between unbiased template and the atlas in commonspace - - `unbiased_to_atlas_inverse_warp/`: inverse of the non-linear transforms for the alignment between unbiased template and the atlas in commonspace - -- `motion_datasink/`: files derivated from motion estimation - - `motion_params_csv/`: contains the 24 motion parameters which can be used as nuisance regressors at the confound correction pipeline stage. - - `FD_csv/`: a CSV file with timescourses for either the mean or maximal framewise displacement (FD) estimations. - - `FD_voxelwise/`: a Nifti image which contains framewise displacement evaluated at each voxel - - `pos_voxelwise/`: a Nifti image which tracks the displacement (derived from the head motion realignment parameters) of each voxel across time - - -## Confound Correction Outputs -Important outputs from confound correction will be found in the `confound_correction_datasink/`: -- `confound_correction_datasink/`: - - `cleaned_timeseries/`: cleaned timeseries after the application of confound correction - - `frame_censoring_mask/`: contains CSV files each recording as a boolean vector which timepoints were censored if frame censoring was applied. - - `aroma_out/`: if `--ica_aroma` is applied, this folder contains outputs from running ICA-AROMA, which includes the MELODIC ICA outputs and the component classification results - - `plot_CR_overfit/`: will contain figures illustrating the variance explained by random regressors during confound correction, and the variance explained by the real regressors after substrating the variance from random regressors. - - -## Analysis Outputs - -Outputs from analyses will be found in the `analysis_datasink/`, whereas outputs relevant to the `--data_diagnosis` are found in `data_diagnosis_datasink/`: -- `analysis_datasink/`: - - `group_ICA_dir/`: complete output from MELODIC ICA, which the melodic_IC.nii.gz Nifti which gives all spatial components, and `report/` folder which includes a HTML visualization. - - `matrix_data_file/`: .pkl file which contains a 2D numpy array representing the whole-brain correlation matrix. If `--ROI_type parcellated` is selected, the row/column indices of the array are matched in increasing order of the atlas ROI label number. - - `matrix_fig/`: .png file which displays the correlation matrix - - `seed_correlation_maps/`: nifti files for seed-based connectivity analysis, where each seed provided in `--seed_list` has an associated voxelwise correlation maps - - `dual_regression_nii/`: the spatial maps from dual regression, which correspond to the linear coefficients from the second regression. The list of 3D spatial maps obtained are concatenated into a 4D Nifti file, where the order of component is consistent with the priors provided in `--prior_maps`. - - `dual_regression_timecourse_csv/`: a CSV file which stores the outputs from the first linear regression during dual regression. This corresponds to a timecourse associated to each prior component from `--prior_maps`. - - `NPR_prior_filename/`: spatial components fitted during NPR - - `NPR_prior_timecourse_csv/`: timecourses associated to each components from NPR_prior_filename - - `NPR_extra_filename/`: the extra spatial components fitted during NPR which were not part of priors - - `NPR_extra_timecourse_csv/`: timecourses associated to each components from NPR_extra_filename -(diagnosis_datasink_target)= -- `data_diagnosis_datasink/`: - - `figure_temporal_diagnosis/`: figure which displays scan-level temporal features from the [spatiotemporal diagnosis](diagnosis_target) - - `figure_spatial_diagnosis/`: figure which displays scan-level spatial features from the [spatiotemporal diagnosis](diagnosis_target) - - `analysis_QC/`: group-level features of data quality from `--data_diagnosis` - - `sample_distributions/`: contains the [distribution plots](dist_plot_target) - - `{analysis}_sample_distribution.png`: the distribution plot for a given network analysis - - `{analysis}_outlier_detection.csv`: a CSV which associates the measures displayed in the distribution plot with corresponding scan IDs - - `parametric_stats/`: [group statistical report](group_stats_target) for analysis quality control (using parametric measures) - - `DR{component #}_QC_maps.png`: The _QC_maps.png files are displaying statistical maps relevant to analysis quality control. The DR refers to dual regression analysis, and the {component #} is relating the file to one of the BOLD components specified in `--prior_bold_idx` - - `DR{component #}_QC_stats.csv`: a follow-up to _QC_maps.png which allows for the quantitative categorization of data quality outcomes as in {cite}`Desrosiers-Gregoire2024-ou` - - `seed_FC{seed #}_QC_maps.png`: same statistical maps as with `DR{component #}_QC_maps.png`, but for seed-based connectivity analysis - - `seed_FC{seed #}_QC_stats.csv`: same measures as with `DR{component #}_QC_maps.png`, but for seed-based connectivity analysis - - `non_parametric_stats/`: same as `parametric_stats/`, but using non-parametric measures - - `temporal_info_csv/`: CSV file containing the data plotted with `figure_temporal_diagnosis/` - - `spatial_VE_nii/`: Nifti file with the confound regression percentage variance explained (R^2) at each voxel - - `CR_prediction_std_nii/`: Nifti file with the confound regression variance explained at each voxel - - `random_CR_std_nii/`: Nifti file with the variance explained from random regressors at each voxel - - `corrected_CR_std_nii/`: Nifti file with the confound regression variance explained at each voxel after removing the variance explained by random regressors - - `temporal_std_nii/`: the standard deviation at each voxel after confound correction - - `GS_cov_nii/`: the covariance of each voxel with the global signal diff --git a/docs/pics/sub-PHG001_ses-3_acq-RARE_T2w_RAS_inho_cor_registration.png b/docs/pics/sub-PHG001_ses-3_acq-RARE_T2w_RAS_inho_cor_registration.png new file mode 100644 index 00000000..6f33bd3a Binary files /dev/null and b/docs/pics/sub-PHG001_ses-3_acq-RARE_T2w_RAS_inho_cor_registration.png differ diff --git a/docs/pics/sub-PHG001_ses-3_acq-RARE_T2w_inho_cor.png b/docs/pics/sub-PHG001_ses-3_acq-RARE_T2w_inho_cor.png new file mode 100644 index 00000000..5697f8d9 Binary files /dev/null and b/docs/pics/sub-PHG001_ses-3_acq-RARE_T2w_inho_cor.png differ diff --git a/docs/pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_inho_cor.png b/docs/pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_inho_cor.png new file mode 100644 index 00000000..2c245938 Binary files /dev/null and b/docs/pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_inho_cor.png differ diff --git a/docs/pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_registration.png b/docs/pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_registration.png new file mode 100644 index 00000000..bcf5917f Binary files /dev/null and b/docs/pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_registration.png differ diff --git a/docs/pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_spatial_diagnosis.png b/docs/pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_spatial_diagnosis.png new file mode 100644 index 00000000..30f10164 Binary files /dev/null and b/docs/pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_spatial_diagnosis.png differ diff --git a/docs/pics/test_dataset_atlas_registration.png b/docs/pics/test_dataset_atlas_registration.png new file mode 100644 index 00000000..dabce302 Binary files /dev/null and b/docs/pics/test_dataset_atlas_registration.png differ diff --git a/docs/preproc_QC.md b/docs/preproc_QC.md deleted file mode 100644 index 70429a55..00000000 --- a/docs/preproc_QC.md +++ /dev/null @@ -1,26 +0,0 @@ -# Preprocessing quality control (QC) - -Several registration operations during preprocessing are prone to fail in accurately aligning images, and it is thus necessary to visually inspect the quality of registration to prevent errors arising from failed alignment, or biased analyses downstream. For this purpose, RABIES generates automatically a set PNG images allowing for efficient visual assessment of key registration steps. These are found in the `{output_folder}/preprocess_QC_report/` folder, which contains several subfolders belonging to different registration step in the pipeline or providing supportive information about the files: - -- `anat_inho_cor/`: intensity inhomogeneities are corrected for prior to important registration operations. This folder allows to assess the quality of the inhomogeneity correction, which is crucial for the performance of downstream registration. The figure is divided in 4 columns, showing 1-the raw image, 2-an initial correction of the image, 3-an overlay of the anatomical mask used to conduct a final correction (by default obtained through a preliminary registration to the commonspace template), and 4-final corrected output. -![](pics/sub-MFC067_ses-1_acq-FLASH_T1w_inho_cor.png) -- `Native2Unbiased/`: alignment between each anatomical image and the generated unbiased template. This registration step controls for the overlap between different scanning sessions. -![](pics/sub-MFC067_ses-1_acq-FLASH_T1w_inho_cor_registration.png) -- `Unbiased2Atlas/`: alignment of the generated unbiased template to the external anatomical template in commonspace. This step ensures proper alignment with the commonspace and the associated brain parcellation. -![](pics/atlas_registration.png) -- `bold_inho_cor/`: same as `anat_inho_cor/`, but conducted on the 3D reference EPI image which is used for estimating the alignment of the EPI. -![](pics/sub-MFC068_ses-1_task-rest_acq-EPI_run-1_bold_inho_cor.png) -- `EPI2Anat/`: shows the alignment of the EPI image to its associated anatomical image of the same scanning session. This step resamples the EPI into native space, and corrects for susceptibility distortions through non-linear registration. An example is shown below: -![](pics/sub-MFC068_ses-1_task-rest_acq-EPI_run-1_bold_registration.png) -- `template_files/`: displays the overlap of the provided external anatomical template with it's associated masks and labels. Allows to validate that proper template files were provided and share those along the RABIES report. -![](pics/template_files.png) -- `temporal_features/`: includes the timecourse of the head motion realignment parameters together with framewise displacement, to observe subject motion. Also includes a spatial map of the signal variability at each voxel and then the temporal signal-to-noise ratio (tSNR). -![](pics/example_temporal_features.png) - - - -```{toctree} ---- -maxdepth: 3 ---- -nested_docs/registration_troubleshoot.md diff --git a/docs/preprocessing.md b/docs/preprocessing.md deleted file mode 100644 index 36fe0873..00000000 --- a/docs/preprocessing.md +++ /dev/null @@ -1,99 +0,0 @@ -# Preprocessing stage -![](pics/preprocessing.png) - -The preprocessing of fMRI scans prior to analysis consists of, at minimum, the anatomical alignment of scans to a common space, head realignment to correct for motion, and the correction of susceptibility distortions arising from the echo-planar imaging (EPI) acquisition of functional scans. The core preprocessing pipeline in RABIES carries each of these steps with state-of-the-art processing tools and techniques. - -To conduct common space alignment, structural images, which were acquired along the EPI scans, are initially corrected for inhomogeneities (**Structural inhomogeneity correction**) and then registered together to allow the alignment of different MRI acquisitions. This registration is conducted by generating an unbiased data-driven template (**Unbiased template generation**) through the iterative non-linear registration of each image to the dataset consensus average, where the average gets updated at each iteration to provide an increasingly representative dataset template (; {cite}`Avants2011-av`). The finalized template after the last iteration provides a representative alignment of each MRI session to a template that shares the acquisition properties of the dataset (e.g. brain shape, FOV, anatomical contrast, ...), making it a stable registration target for cross-subject alignment. This newly-generated unbiased template is then itself registered to an external reference atlas to provide both an anatomical segmentation and a common space comparable across studies defined from the provided reference atlas (**Atlas registration**). - -The remaining preprocessing involves the EPI image. A volumetric EPI image is first derived using a trimmed mean across the EPI frames, after an initial motion realignment step (**3D EPI generation**). Using this volumetric EPI as a target, the head motion parameters are estimated by realigning each EPI frame to the target using a rigid registration (**Head motion estimation**). To correct for EPI susceptibility distortions, the volumetric EPI is first subjected to an inhomogeneity correction step (**Functional inhomogeneity correction**), and then registered non-linearly to the anatomical scan from the same MRI session, which allows to calculate the required geometrical transforms for recovering brain anatomy {cite}`Wang2017-ci` (**Susceptibility distortion estimation**). Finally, after calculating the transformations required to correct for head motion and susceptibility distortions, both transforms are concatenated into a single resampling operation (avoiding multiple resampling) which is applied at each EPI frame, generating the preprocessed EPI timeseries in native space {cite}`Esteban2019-rs` (**Frame-wise resampling**). Preprocessed timeseries in common space are also generated by further concatenating the transforms allowing resampling to the reference atlas. - -The workflow of the RABIES preprocessing pipeline is summarized in the diagram above, and each preprocessing module is further described below. - -## Structural inhomogeneity correction -![](pics/sub-MFC067_ses-1_acq-FLASH_T1w_inho_cor.png) -**Figure:** displays steps of inhomogeneity correction for the structural image. - -### rabies.preprocess_pkg.inho_correction.init_inho_correction_wf [[source code](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/inho_correction.py)] - -```{literalinclude} ../rabies/preprocess_pkg/inho_correction.py -:start-after: inho_correction_head_start -:end-before: inho_correction_head_end -``` - -## Common space alignment (i.e. Unbiased template generation + Atlas registration) -![](pics/sub-MFC067_ses-1_acq-FLASH_T1w_inho_cor_registration.png) -**Figure:** displays the overlap between a structural scan (top) and the dataset-generated unbiased template (bottom). -![](pics/atlas_registration.png) -**Figure:** displays the overlap between the unbiased template (top) and the reference atlas template (bottom). - -### rabies.preprocess_pkg.commonspace_reg.init_commonspace_reg_wf [[source code](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/commonspace_reg.py)] -```{literalinclude} ../rabies/preprocess_pkg/commonspace_reg.py -:start-after: commonspace_wf_head_start -:end-before: commonspace_wf_head_end -``` - -(3D_EPI_target)= -## 3D EPI generation -### rabies.preprocess_pkg.bold_ref.init_bold_reference_wf [[source code](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/bold_ref.py)] - -```{literalinclude} ../rabies/preprocess_pkg/bold_ref.py -:start-after: gen_bold_ref_head_start -:end-before: gen_bold_ref_head_end -``` - -## Head motion estimation -![](pics/example_motion_parameters.png) -**Figure:** example of the 6 motion parameters. - -### rabies.preprocess_pkg.hmc.init_bold_hmc_wf [[source code](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/hmc.py)] - -```{literalinclude} ../rabies/preprocess_pkg/hmc.py -:start-after: hmc_wf_head_start -:end-before: hmc_wf_head_end -``` - -### rabies.preprocess_pkg.hmc.EstimateMotionParams [[source code](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/hmc.py)] - -```{literalinclude} ../rabies/preprocess_pkg/hmc.py -:start-after: motion_param_head_start -:end-before: motion_param_head_end -``` - -## Functional inhomogeneity correction -![](pics/sub-MFC068_ses-1_task-rest_acq-EPI_run-1_bold_inho_cor.png) -**Figure:** displays steps of inhomogeneity correction for the volumetric EPI. - -The workflow is the same as the **structural inhomogeneity correction**. - -## Susceptibility distortion estimation -![](pics/sub-MFC068_ses-1_task-rest_acq-EPI_run-1_bold_registration.png) -**Figure:** displays the overlap between the volumetric EPI (top) and structural image (bottom). - -### rabies.preprocess_pkg.registration.init_cross_modal_reg_wf [[source code](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/registration.py)] - -```{literalinclude} ../rabies/preprocess_pkg/registration.py -:start-after: cross_modal_reg_head_start -:end-before: cross_modal_reg_head_end -``` - -## Frame-wise resampling - -### rabies.preprocess_pkg.resampling.init_bold_preproc_trans_wf [[source code](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/resampling.py)] - -```{literalinclude} ../rabies/preprocess_pkg/resampling.py -:start-after: bold_resampling_head_start -:end-before: bold_resampling_head_end -``` - -## Resampling of masks and labels - -### rabies.preprocess_pkg.resampling.init_mask_preproc_trans_wf [[source code](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/resampling.py)] - -```{literalinclude} ../rabies/preprocess_pkg/resampling.py -:start-after: mask_resampling_head_start -:end-before: mask_resampling_head_end -``` - -## Adapted workflow without structural scans (i.e. --bold_only) -Structural scans are recommended, but not required to complete preprocessing with RABIES. An alternative workflow is also implemented to preprocess a input dataset which contains only EPI functional images, and can be selected with the `--bold_only` option. In this alternative workflow, the volumetric EPI corrected for inhomogeneity during **Functional inhomogeneity correction** replaces the structural image for the purpose of common space alignment, and is thus used for generating the unbiased template, in turn, this template is registered to the reference atlas. This final registration to the atlas accounts for estimation of susceptibility distortions instead of the registration to a structural image from the same MRI session. -If using the RABIES default mouse atlas, the default template is changed to a EPI reference template, which offers a more robust target for EPI registration than a structural image as reference template. diff --git a/docs/reference/bibliography.md b/docs/reference/bibliography.md new file mode 100644 index 00000000..7fb4808b --- /dev/null +++ b/docs/reference/bibliography.md @@ -0,0 +1,6 @@ +# Bibliography + +Every work cited in this documentation. + +```{bibliography} ../_static/refs.bib +``` diff --git a/docs/reference/bids_inputs.md b/docs/reference/bids_inputs.md new file mode 100644 index 00000000..94f7f384 --- /dev/null +++ b/docs/reference/bids_inputs.md @@ -0,0 +1,94 @@ +# Input data requirements + +The input dataset must be organised according to the +[BIDS data structure](https://bids.neuroimaging.io/) {cite}`Gorgolewski2016-zm`. + +## How RABIES traverses a dataset + +RABIES iterates through every subject found to contain a functional file, and +within each subject through the sessions and runs present. + +When anatomical scans are used — that is, when `--bold_only` is not set — each +functional scan is matched to one anatomical scan **from the same subject and +session**. + +## Scan identification + +```{list-table} +:header-rows: 1 +:widths: 30 30 40 + +* - Image type + - Default BIDS suffixes + - Parameter +* - Functional + - `bold`, `cbv` + - `--bids_filter`, key `func` +* - Structural + - `T1w`, `T2w` + - `--bids_filter`, key `anat` +``` + +Files matching neither filter are ignored. The default filter is equivalent to: + +```{code-block} json +:caption: Default value of --bids_filter + +{ + "func": { + "suffix": ["bold", "cbv"] + }, + "anat": { + "suffix": ["T1w", "T2w"] + } +} +``` + +```{seealso} +[How to select which scans get processed](../how_to/select_scans.md) for +customising the filter and for selecting individual scans. +``` + +## Image orientation + +RABIES expects images in the NIfTI standard RAS+ orientation +(Right–Anterior–Superior). Incorrectly oriented images are a common source of +registration failures — see +[How to check image orientation](../how_to/check_orientation.md). + +## Example dataset + +The [RABIES example dataset](http://doi.org/10.5281/zenodo.8349029) +(`test_dataset.zip`) has the following structure: + +```{code-block} text +:caption: Two subjects, one session each, with paired anatomical and functional scans + +test_dataset/ +├── sub-PHG001 +│ └── ses-3 +│ ├── anat +│ │ ├── sub-PHG001_ses-3_acq-RARE_T2w.json +│ │ └── sub-PHG001_ses-3_acq-RARE_T2w.nii.gz +│ └── func +│ ├── sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold.json +│ └── sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold.nii.gz +└── sub-PHG002 + └── ses-3 + ├── anat + │ ├── sub-PHG002_ses-3_acq-RARE_T2w.json + │ └── sub-PHG002_ses-3_acq-RARE_T2w.nii.gz + └── func + ├── sub-PHG002_ses-3_task-rest_acq-EPI_run-1_bold.json + └── sub-PHG002_ses-3_task-rest_acq-EPI_run-1_bold.nii.gz + +8 directories, 8 files +``` + +This is the dataset used in [the tutorial](../tutorials/first_run.md). + +## Format conversion + +Conversion from Bruker raw format to NIfTI can be handled with +[BrkRaw](https://brkraw.github.io/). The CoBrALab maintains +[notes on the conversion](https://github.com/CoBrALab/documentation/wiki/bruker2nifti-conversion). diff --git a/docs/reference/cli.md b/docs/reference/cli.md new file mode 100644 index 00000000..1146dacf --- /dev/null +++ b/docs/reference/cli.md @@ -0,0 +1,102 @@ +# Command line interface + +RABIES runs from a terminal. The software is divided into three processing +stages, and one of them must be selected on every invocation: + +```text +rabies [execution options] {preprocess,confound_correction,analysis} [stage options] input output +``` + +```{list-table} +:header-rows: 1 +:widths: 25 25 50 + +* - Stage + - Takes as input + - Produces +* - `preprocess` + - A BIDS directory + - Motion-corrected, distortion-corrected timeseries aligned to commonspace +* - `confound_correction` + - A `preprocess` output directory + - Cleaned timeseries +* - `analysis` + - A `confound_correction` output directory + - Connectivity estimates and data quality reports +``` + +Execution options — parallel processing, memory management, scan selection — +are given **before** the stage name. Stage options are given after it. + +```{note} +The `--help` output below is generated at documentation build time from the +RABIES version these docs were built for. Run `rabies --help` locally to see +the options for the version you have installed. +``` + +## Argument syntax + +RABIES uses diverse argument shapes: + +Flags +: `--apply_STC`. Present or absent, no value. + +Single values +: `--TR 1.2`. The value follows the parameter name after a space. + +Lists +: `--nuisance_regressors WM_signal CSF_signal mot_6`. Values follow the + parameter name, separated by spaces. + +Key-value groups +: `--commonspace_reg masking=true,template_registration=SyN`. One parameter + takes several settings as comma-separated `key=value` pairs, with no spaces. + +## `rabies --help` + +Execution options that apply to every stage: parallel processing plugin, thread +and memory limits, scan inclusion and exclusion, output data type and +interpolation. + +:::{dropdown} rabies --help +:icon: terminal + +```{program-output} rabies --help +``` +::: + +## `rabies preprocess --help` + +Input selection, image corrections, registration options, resampling, slice +timing correction, and the reference atlas files. + +:::{dropdown} rabies preprocess --help +:icon: terminal + +```{program-output} rabies preprocess --help +``` +::: + +## `rabies confound_correction --help` + +Frame censoring, detrending, ICA-AROMA, frequency filtering, confound +regression, intensity scaling and smoothing. + +:::{dropdown} rabies confound_correction --help +:icon: terminal + +```{program-output} rabies confound_correction --help +``` +::: + +## `rabies analysis --help` + +Seed-based connectivity, whole-brain connectivity matrices, group ICA, dual +regression, and the `--data_diagnosis` quality reports. + +:::{dropdown} rabies analysis --help +:icon: terminal + +```{program-output} rabies analysis --help +``` +::: diff --git a/docs/reference/index.md b/docs/reference/index.md new file mode 100644 index 00000000..f0b2b700 --- /dev/null +++ b/docs/reference/index.md @@ -0,0 +1,37 @@ +# Reference + +Descriptive, factual material to consult while you work. Reference pages +describe how RABIES behaves; they do not tell you what to do. For that, see the +[how-to guides](../how_to/index.md). + +## Running RABIES + +```{toctree} +--- +maxdepth: 1 +--- +bids_inputs +template_files +cli +``` + +## What RABIES produces + +```{toctree} +--- +maxdepth: 1 +--- +outputs +qc_outputs +metrics +``` + +## Internals + +```{toctree} +--- +maxdepth: 1 +--- +workflows +bibliography +``` diff --git a/docs/reference/metrics.md b/docs/reference/metrics.md new file mode 100644 index 00000000..c1bea3f3 --- /dev/null +++ b/docs/reference/metrics.md @@ -0,0 +1,240 @@ +(metrics_target)= + +# Metric definitions + +Precise definitions of every quantity RABIES computes and reports. + +```{seealso} +For where these values appear, see [Output files](outputs.md). For how to +interpret them, see [Data quality assessment](../explanation/data_quality.md). +``` + +(regressor_target)= + +## Nuisance regressors + +Selected with `--nuisance_regressors` at the confound correction stage. + +**mot_6** +: 3 rotations (Euler angles in radians) and 3 translations (in mm) measured from + the rigid-body head realignment algorithm. + +**mot_24** +: The 6 motion parameters together with their temporal derivatives, plus 12 + additional parameters obtained by taking the squared terms — the Friston 24 + parameters {cite}`Friston1996-sa`: + + $$ + mot24_t = [mot6_t,(mot6_t-mot6_{t-1}),(mot6_t)^2,(mot6_t-mot6_{t-1})^2] + $$ + + with $mot24_t$ representing the list of 24 regressors for timepoint $t$. + +**WM/CSF/vascular/global signal** +: The mean signal computed within the corresponding brain mask (WM, CSF, + vascular or whole-brain). + +**aCompCor_percent** +: Principal component timecourses derived from timeseries within the combined + WM and CSF masks — the aCompCor technique {cite}`Muschelli2014-vi`. From the + timeseries within the WM/CSF masks $Y_{WM/CSF}$, a principal component + analysis (PCA) decomposition is conducted to derive + + $$ + Y_{WM/CSF} = W_{aCompCor}C^T + $$ + + with $C$ a set of spatial principal components and $W$ their associated + loadings across time. The first components explaining 50% of the variance are + kept, and their loadings $W_{aCompCor}$ provide the aCompCor nuisance + regressors. The PCA is conducted on partially cleaned timeseries after + step 4 of the confound correction pipeline (see [The confound correction workflow](../explanation/confound_correction.md)). + +**aCompCor_5** +: As **aCompCor_percent**, but the first 5 components are kept instead of a set + explaining 50% of the variance. + +## Temporal scan diagnosis + +(mot6_target)= + +**Head motion translation and rotation parameters** +: 3 rotations (Euler angles in radians) and 3 translations (in mm) measured from + the rigid-body head realignment algorithm. + +(FD_target)= + +**Framewise displacement** +: For each timepoint, this is the displacement between + the current and the previous frame, averaged over all brain voxels. + For each brain voxel within the referential + space for head realignment (i.e. the [3D EPI](wf_bold_ref) provided as reference + for realignment) and for each timepoint, the inverse transform of the head + motion parameters from the corresponding timepoint is applied to obtain the + voxel position pre-motion correction. Framewise displacement is then computed + for each voxel as the Euclidean distance between the pre-motion-correction + positions for the current and previous timepoints. The mean framewise + displacement $FD_t$ at timepoint $t$ is therefore + + $$ + FD_t = \frac{1}{n}\sum_{i=1}^{n}\sqrt{(x_{i,t}-x_{i,t-1})^2+(y_{i,t}-y_{i,t-1})^2+(z_{i,t}-z_{i,t-1})^2} + $$ + + using the 3D $x$, $y$ and $z$ spatial coordinates in mm for timepoints $t$ + and $t-1$ and voxel indices $i$. Framewise displacement for the first frame, + which has no past timepoint, is set to 0. + +(DVARS_target)= + +**DVARS** +: The estimation of temporal shifts in global signal at each timepoint, + measured as the root-mean-square of the timeseries' temporal derivative + + $$ + DVARS_t = \sqrt{\frac{1}{n}\sum_{i=1}^{n}(Y_{i,t}-Y_{i,t-1})^2} + $$ + + where $Y_{i,t}$ is the BOLD signal in brain voxel $i$ at timepoint $t$. The + first timepoint is set to 0, having no previous timepoint. + +**Framewise distance from mean** +: The mean square error (MSE) between each frame and the average EPI, with the + average computed as the tri-mean across time voxelwise. The input is the + resampled timeseries output from the `preprocess` stage. + +**Whole-brain/Edge/WM/CSF mask** +: The mean signal across a given brain mask: whole-brain (the global signal), + WM, CSF or brain edge. + +**CRvar** +: The variance estimated and removed by nuisance regression at each timepoint, computed as + the root-mean square (RMS) of $Y_{CR}$ across voxels, where $Y_{CR}$ is the + [predicted confound timeseries](CR_target). + +**CR R2** +: The proportion of variance explained, and removed, by nuisance regression. + Obtained with $CR_{R^2}= 1-\frac{Var(\hat{Y_t})}{Var(Y_t)}$ for a given frame at timepoint $t$, + where $Y_t$ and $\hat{Y_t}$ are the BOLD frame at timepoint $t$ pre- and post-regression, and + $Var(x) = \frac{1}{n}\sum_{i=1}^{n}(x_i - \mu_x)^2$ is the variance with + $\mu$ the mean. + +**Mean amplitude of network VS confound sources** +: A set of timecourses averaged as $\frac{1}{n}\sum_{i=1}^{n}|X_i|$, where + $X_i$ is timecourse $i$. The timecourses correspond to one of: + + - *DR confounds*: timecourses from the first stage of dual regression, using + the confound components provided to `--prior_confound_idx` + - *DR networks*: network timecourses from the first stage of dual regression + as specified with `--prior_bold_idx` + - *SBC networks*: network timecourses derived from the seeds provided in + `--seed_list` + +## Spatial scan diagnosis + +**BOLDSD** +: The temporal standard deviation computed for each voxel from the cleaned BOLD + timeseries. + +**CRSD** +: The temporal standard deviation computed for each voxel from the predicted + confound timeseries during confound regression, i.e. [$Y_{CR}$](CR_target). + +**CR R2** +: The proportion of variance explained by confound regression at each voxel. + Obtained with $CR_{R^2}= 1-\frac{Var(\hat{Y_v})}{Var(Y_v)}$ for a given voxel $v$, where + $Y_v$ and $\hat{Y_v}$ are the timeseries pre- and post-regression for that voxel, and + $Var(x) = \frac{1}{n}\sum_{i=1}^{n}(x_i - \mu_x)^2$ is the variance of $x$ + with $\mu$ the mean. + +**Global signal covariance (GScov)** +: The covariance between the mean whole-brain signal (i.e. global signal) timecourse + and the timeseries at each voxel. + +**DR network X** +: The linear coefficients resulting from the + [second regression with dual regression](DR_target), corresponding to a + network amplitude map, for the Xth network specified with `--prior_bold_idx`. + +**SBC network X** +: The voxelwise correlation coefficients (Pearson's r) estimated with + seed-based connectivity, for the Xth seed provided in `--seed_list`. + +(dist_plot_metrics)= + +## Distribution plot + +**Network amplitude** +: The overall network amplitude, summarised by computing the L2-norm across a + network connectivity map from a subject-level analysis. Such a map can be + derived from seed-based correlation, or correspond to the linear coefficients + from the [second regression ${\beta}_{SM}$](DR_target) for dual regression. + +**Network specificity** +: The network map for a given analysis (seed-based or dual regression) and a given + fMRI scan is compared relative to a reference network map to assess specificity. + The reference map is either the original ICA component for dual regression, a + manually provided input map with `--seed_prior_list` for seed-based analysis, + or the dataset average itself if selecting `--group_avg_prior`. + To compute specificity, both the reference and individual scan network maps + are thresholded to retain the top X% of voxels with + highest connectivity, X% being defined by `--brainmap_percent_threshold`, and + the overlap of the thresholded area is computed using Dice overlap. + +```{important} +If the reference network map was generated from the dataset average, it is +important to validate that the average can indeed provide an adequate +representation of the expected network connectivity. If the reference does +not represent a network, then the network specificity metric is meaningless. +The average network used for those computations can be visualised in the +[group statistical report](group_stats_target). +``` + +**Dual regression confound correlation** +: The timecourse for a single network, from a seed or from dual regression, is + correlated with the timecourse from each confound component (provided using + `--prior_confound_idx`) modelled through dual regression. The absolute mean + correlation is then computed to obtain the average amplitude of confound + correlations for that network analysis. + +**FD-DVARS corr.** +: For each scan, the correlation between the framewise displacement timecourse + and the DVARS timecourse **computed post-confound correction** — this is not + the DVARS plotted in the temporal diagnosis figure. Censored timeframes are + excluded from both timecourses, and DVARS is recomputed after applying + confound correction, so this metric represents *residual* associations + between spontaneous motion and the cleaned global signal fluctuations. + +**Total $CR_{SD}$** +: The total standard deviation across the + [predicted confound timeseries $Y_{CR}$](CR_target). + +**Mean framewise displacement** +: The mean framewise displacement computed across time, including only frames + remaining after the censoring applied for confound correction. + +**Temporal degrees of freedom** +: The degrees of freedom remaining after confound correction: + + ```text + tDOF = Original number of timepoints + - Number of censored timepoints + - Number of AROMA components removed + - Number of nuisance regressors + ``` + +(group_QC_metrics)= + +## Group statistical QC report + +These quantities are stored inside the `analysis_QC/*_stats/*_QC_stats.csv` output files. + +**Specificity of network variability** +: As with network specificity in the distribution plot, the network variability + map and the corresponding reference network map are thresholded to include + the top X% of voxels (X% defined by `--brainmap_percent_threshold`), and the + overlap is estimated using Dice overlap. + +**Mean confound correlation** +: For each confound correlation map ($CR_{SD}$, mean FD or tDOF), the mean is + computed across voxels within the thresholded area of the reference network + map, giving a mean correlation within the network's core region. diff --git a/docs/reference/outputs.md b/docs/reference/outputs.md new file mode 100644 index 00000000..5ed90d63 --- /dev/null +++ b/docs/reference/outputs.md @@ -0,0 +1,180 @@ +# Output files + +Every processing stage writes its important outputs into `datasink/` folders, +created inside the output directory given at execution. This page describes +every file produced. + +```{note} +The output directory also contains the Nipype working directory and the +serialised workflow state (`.pkl`). Those are internal, but they are what the +next stage reads to locate the previous stage's files — do not move or delete +them between stages. +``` + +## Preprocessing outputs + +Preprocessing writes five datasinks: `anat_datasink/`, `bold_datasink/`, +`unbiased_template_datasink/`, `transforms_datasink/` and `motion_datasink/`. +The quality control images are written separately to `preprocess_QC_report/`, +described in [Preprocessing QC outputs](qc_outputs.md). + +### `anat_datasink/` + +Inhomogeneity-corrected anatomical scans. + +- `anat_preproc/`: anatomical scans after inhomogeneity correction + +### `bold_datasink/` + +All outputs related to the functional scans. Files are resampled either onto +the native space or the commonspace of the EPI. Native space outputs are +resampled over the anatomical scan from the corresponding MRI session; +commonspace outputs are resampled over the reference atlas. The original EPI +voxel resolution is unchanged during resampling unless specified otherwise in +the RABIES command. + +**Native space** + +- `native_bold/`: preprocessed EPI timeseries resampled to native space +- `native_brain_mask/`: brain mask in native space +- `native_WM_mask/`: WM mask in native space +- `native_CSF_mask/`: CSF mask in native space +- `native_labels/`: atlas labels in native space +- `native_bold_ref/`: a volumetric 3D EPI average generated from the 4D `native_bold/` + +**Commonspace** + +- `commonspace_bold/`: preprocessed EPI timeseries resampled to commonspace +- `commonspace_mask/`: brain mask in commonspace +- `commonspace_WM_mask/`: WM mask in commonspace +- `commonspace_CSF_mask/`: CSF mask in commonspace +- `commonspace_vascular_mask/`: vascular mask in commonspace +- `commonspace_labels/`: atlas labels in commonspace +- `commonspace_resampled_template/`: the commonspace anatomical template, resampled to the EPI's dimensions + +**Inputs and intermediates** + +- `input_bold/`: the raw EPI scans provided as inputs in the BIDS data folder +- `initial_bold_ref/`: the initial volumetric 3D EPI average generated from the 4D `input_bold/` +- `raw_brain_mask/`: brain mask resampled onto the 4D `input_bold/` +- `inho_cor_bold/`: the volumetric 3D EPI (`initial_bold_ref/`) after inhomogeneity correction, later used for registration of the EPI +- `inho_cor_bold_warped2anat/`: `inho_cor_bold` after co-registration to the associated anatomical image (`anat_preproc/`) +- `std_map_preprocess/`: the temporal standard deviation at each voxel of `commonspace_bold/` +- `tSNR_map_preprocess/`: the temporal signal-to-noise ratio (tSNR) of `commonspace_bold/` + +### `unbiased_template_datasink/` + +Outputs from the generation of the unbiased template using +[optimized_antsMultivariateTemplateConstruction](https://github.com/CoBrALab/optimized_antsMultivariateTemplateConstruction). +The unbiased template is the average of all anatomical scans (or functional +scans with `--bold_only`) after their alignment. + +- `unbiased_template/`: the unbiased template generated from the input dataset scans +- `warped_unbiased_template/`: the unbiased template, registered to the reference atlas in commonspace + +### `transforms_datasink/` + +All transform files for resampling between spaces. + +The `bold_to_anat` registration transforms the raw EPI to overlap with the +anatomical image, correcting susceptibility distortions, which defines native +space. The `native_to_unbiased` registration overlaps every scan onto the +generated unbiased template. The `unbiased_to_atlas` registration aligns the +unbiased template with the reference atlas, which defines commonspace. + +- `bold_to_anat_affine/`: affine transforms from the EPI co-registration to the anatomical image +- `bold_to_anat_warp/`: non-linear transforms from the EPI co-registration to the anatomical image +- `bold_to_anat_inverse_warp/`: inverse of `bold_to_anat_warp/` +- `native_to_unbiased_affine/`: affine transforms for the alignment between native space and the unbiased template +- `native_to_unbiased_warp/`: non-linear transforms for the same alignment +- `native_to_unbiased_inverse_warp/`: inverse of `native_to_unbiased_warp/` +- `unbiased_to_atlas_affine/`: affine transforms for the alignment between the unbiased template and the atlas in commonspace +- `unbiased_to_atlas_warp/`: non-linear transforms for the same alignment +- `unbiased_to_atlas_inverse_warp/`: inverse of `unbiased_to_atlas_warp/` + +### `motion_datasink/` + +Files derived from motion estimation. + +- `motion_params_csv/`: the 24 motion parameters, usable as nuisance regressors at the confound correction stage +- `FD_csv/`: a CSV with timecourses for either the mean or maximal [framewise displacement](FD_target) estimations +- `FD_voxelwise/`: a NIfTI image containing framewise displacement evaluated at each voxel +- `pos_voxelwise/`: a NIfTI image tracking the displacement of each voxel across time, derived from the head motion realignment parameters + +## Confound correction outputs + +### `confound_correction_datasink/` + +- `cleaned_timeseries/`: cleaned timeseries after the application of confound correction +- `frame_censoring_mask/`: CSV files recording, as a boolean vector, which timepoints were censored, if frame censoring was applied +- `aroma_out/`: outputs from running ICA-AROMA if `--ica_aroma` is applied, including the MELODIC ICA outputs and the component classification results +- `plot_CR_overfit/`: figures illustrating the variance explained by random regressors during confound correction, and the variance explained by the real regressors after subtracting the variance from random regressors + +## Analysis outputs + +### `commonspace_analysis_datasink/` and `nativespace_analysis_datasink/` + +Which of these appears depends on the space the cleaned timeseries were +produced in: + +```{list-table} +:header-rows: 1 +:widths: 45 55 + +* - Datasink + - Present when +* - `commonspace_analysis_datasink/` + - confound correction produced commonspace timeseries, or `--resample_to_commonspace` was passed, or `--group_ica apply=true` was used +* - `nativespace_analysis_datasink/` + - confound correction produced nativespace timeseries +``` + +Both contain the same set of files, listed below. `group_ICA_dir/` is always +written to `commonspace_analysis_datasink/`, since group ICA requires +commonspace alignment. + +- `group_ICA_dir/`: complete output from MELODIC ICA, including the `melodic_IC.nii.gz` NIfTI giving all spatial components and a `report/` folder with an HTML visualisation +- `matrix_data_file/`: a `.pkl` file containing a 2D NumPy array representing the whole-brain correlation matrix. With `--ROI_type parcellated`, the row and column indices are matched in increasing order of the atlas ROI label number +- `matrix_fig/`: a `.png` displaying the correlation matrix +- `seed_correlation_maps/`: NIfTI files for [seed-based connectivity](SBC_target), one voxelwise correlation map per seed provided in `--seed_list` +- `dual_regression_nii/`: the spatial maps from [dual regression](DR_target), corresponding to the linear coefficients from the second regression. The 3D spatial maps are concatenated into a 4D NIfTI, with component order consistent with the priors provided in `--prior_maps` +- `dual_regression_timecourse_csv/`: a CSV storing the outputs from the first linear regression during dual regression — one timecourse per prior component from `--prior_maps` +- `NPR_prior_filename/`: spatial components fitted during NPR +- `NPR_prior_timecourse_csv/`: timecourses associated with each component in `NPR_prior_filename/` +- `NPR_extra_filename/`: the extra spatial components fitted during NPR which were not part of the priors +- `NPR_extra_timecourse_csv/`: timecourses associated with each component in `NPR_extra_filename/` + +(diagnosis_datasink_target)= + +### `data_diagnosis_datasink/` + +Produced when `--data_diagnosis` is selected. See +[How to assess data quality](../how_to/assess_data_quality.md) for how to use +these. + +- `figure_temporal_diagnosis/`: scan-level temporal features from the [spatiotemporal diagnosis](diagnosis_target) +- `figure_spatial_diagnosis/`: scan-level spatial features from the [spatiotemporal diagnosis](diagnosis_target) +- `temporal_info_csv/`: CSV containing the data plotted in `figure_temporal_diagnosis/` +- `spatial_VE_nii/`: NIfTI with the confound regression percentage variance explained ($R^2$) at each voxel +- `CR_prediction_std_nii/`: NIfTI with the confound regression variance explained at each voxel +- `random_CR_std_nii/`: NIfTI with the variance explained from random regressors at each voxel +- `corrected_CR_std_nii/`: NIfTI with the confound regression variance explained at each voxel after removing the variance explained by random regressors +- `temporal_std_nii/`: the standard deviation at each voxel after confound correction +- `GS_cov_nii/`: the covariance of each voxel with the global signal + +`analysis_QC/` holds the group-level features of data quality: + +- `sample_distributions/`: the [distribution plots](dist_plot_target) + - `{analysis}_sample_distribution.png`: the distribution plot for a given network analysis + - `{analysis}_outlier_detection.csv`: a CSV associating the measures displayed in the distribution plot with the corresponding scan IDs +- `parametric_stats/`: the [group statistical report](group_stats_target) for analysis quality control, using parametric measures + - `DR{component #}_QC_maps.png`: statistical maps relevant to analysis quality control. `DR` refers to dual regression analysis, and `{component #}` relates the file to one of the BOLD components specified in `--prior_bold_idx` + - `DR{component #}_QC_stats.csv`: a follow-up to `_QC_maps.png` allowing quantitative categorisation of data quality outcomes as in {cite}`Desrosiers-Gregoire2024-ou` + - `seed_FC{seed #}_QC_maps.png`: the same statistical maps, for seed-based connectivity analysis + - `seed_FC{seed #}_QC_stats.csv`: the same measures, for seed-based connectivity analysis +- `non_parametric_stats/`: as `parametric_stats/`, but using non-parametric measures + +```{seealso} +[Metric definitions](metrics.md) gives the precise computation behind every +quantity named on this page. +``` diff --git a/docs/reference/qc_outputs.md b/docs/reference/qc_outputs.md new file mode 100644 index 00000000..4d4b7fb2 --- /dev/null +++ b/docs/reference/qc_outputs.md @@ -0,0 +1,160 @@ +# Preprocessing QC outputs + +Several registration operations during preprocessing can fail to align images +accurately. RABIES automatically generates a set of PNG images allowing +efficient visual assessment of the key registration steps, so that failed +alignment is caught before it biases downstream analyses. + +The images are written to `{output_folder}/preprocess_QC_report/`, one +subfolder per registration step or per supporting file group. + +```{important} +Inspect this report on every run. A failed registration does not raise an +error — it produces plausible-looking output that is wrong. +``` + +## Which folders you get + +Some folder names carry the name of the workflow that produced them as a +prefix, and which commonspace folders exist depends on `--bold_only`, +`--commonspace_reg`, `--anat_robust_inho_cor` and `--bold_robust_inho_cor`. + +```{list-table} +:header-rows: 1 +:widths: 45 55 + +* - Folder + - Present when +* - `anat_inho_cor/` + - structural scans are used, i.e. not `--bold_only` +* - `bold_inho_cor/` + - always +* - `commonspace_reg_wf.Anat2Unbiased/` + - `fast_commonspace=false` (the default) +* - `commonspace_reg_wf.Unbiased2Atlas/` + - `fast_commonspace=false` (the default) +* - `commonspace_reg_wf.unbiased_template_masking/` + - `fast_commonspace=false` and `masking=true` +* - `commonspace_reg_wf.Anat2Atlas/` + - `fast_commonspace=true` +* - `anat_robust_inho_cor_template.*/` + - `--anat_robust_inho_cor apply=true` +* - `bold_robust_inho_cor_template.*/` + - `--bold_robust_inho_cor apply=true` +* - `EPI2Anat/` + - structural scans are used, i.e. not `--bold_only` +* - `template_files/` + - always +* - `temporal_features/` + - always +``` + +The `*_robust_inho_cor_template.*` folders contain the same registration +figures as the `commonspace_reg_wf.*` ones, for the template built during the +robust inhomogeneity correction pass rather than for the main commonspace +registration. + +## `anat_inho_cor/` + +Quality of the intensity inhomogeneity correction applied to the structural +image, which is performed before the important registration operations and is +crucial for their performance. + +Each figure has 4 columns: **1** the raw image, **2** an initial correction of +the image, **3** an overlay of the anatomical mask used for the final +correction (by default obtained through a preliminary registration to the +commonspace template), and **4** the final corrected output. + +```{figure} ../pics/sub-MFC067_ses-1_acq-FLASH_T1w_inho_cor.png +:alt: Four-column figure showing the stages of structural inhomogeneity correction + +Structural inhomogeneity correction. +``` + +## `bold_inho_cor/` + +The same as `anat_inho_cor/`, but for the 3D reference EPI image used to +estimate the alignment of the EPI. + +```{figure} ../pics/sub-MFC068_ses-1_task-rest_acq-EPI_run-1_bold_inho_cor.png +:alt: Four-column figure showing the stages of EPI inhomogeneity correction + +Functional inhomogeneity correction. +``` + +## `commonspace_reg_wf.Anat2Unbiased/` + +Alignment between each anatomical image and the generated unbiased template. +This registration controls the overlap between different scanning sessions. + +```{figure} ../pics/sub-MFC067_ses-1_acq-FLASH_T1w_inho_cor_registration.png +:alt: Overlap between a structural scan and the dataset-generated unbiased template + +Structural scan (top) against the unbiased template (bottom). The red lining is +an automated edge-detection computed from the **opposite** image - the lining on +the top image is computed from the bottom one and vice-versa. This enables +conveniently comparing edge-alignment. +``` + +## `commonspace_reg_wf.Unbiased2Atlas/` + +Alignment of the generated unbiased template to the external anatomical +template in commonspace. This ensures proper alignment with the commonspace and +its associated brain parcellation. + +```{figure} ../pics/atlas_registration.png +:alt: Overlap between the unbiased template and the reference atlas template + +Unbiased template (top) against the reference atlas template (bottom). The red +lining on each image is the automated edge-detection computed from the other image. +``` + +## `commonspace_reg_wf.Anat2Atlas/` + +Produced instead of the two folders above when `fast_commonspace=true`, which +skips the unbiased template and registers each scan directly to the reference +atlas. It shows that direct alignment, one figure per scan. + +## `EPI2Anat/` + +Alignment of the EPI image to the anatomical image from the same scanning +session. This step resamples the EPI into native space and corrects +susceptibility distortions through non-linear registration. + +```{figure} ../pics/sub-MFC068_ses-1_task-rest_acq-EPI_run-1_bold_registration.png +:alt: Overlap between the volumetric EPI and the structural image + +Volumetric EPI (top) against the structural image (bottom). The red lining on each +image is the automated edge-detection computed from the other image. +``` + +## `template_files/` + +Overlap of the provided external anatomical template with its associated masks +and labels. Use it to confirm the correct template files were provided, and +share it alongside the RABIES report. + +```{figure} ../pics/template_files.png +:alt: The anatomical template overlaid with its masks and labels + +Template, masks and labels. +``` + +## `temporal_features/` + +The timecourse of the head motion realignment parameters together with +framewise displacement, showing subject motion. Also includes a spatial map of +the signal variability at each voxel, and the temporal signal-to-noise ratio +(tSNR). + +```{figure} ../pics/example_temporal_features.png +:alt: Motion parameter timecourses, framewise displacement, variability and tSNR maps + +Temporal features. +``` + +```{seealso} +- [How to troubleshoot registration](../how_to/troubleshoot_registration.md) — what to change when these figures show a failure +- [Metric definitions](metrics.md) — how framewise displacement and the other quantities are computed +- [Data quality assessment](../explanation/data_quality.md) — the separate, analysis-stage quality reports +``` diff --git a/docs/reference/template_files.md b/docs/reference/template_files.md new file mode 100644 index 00000000..29bbe534 --- /dev/null +++ b/docs/reference/template_files.md @@ -0,0 +1,25 @@ +# Built-in commonspace template and atlas + +## The DSURQE mouse atlas + +By default, RABIES uses the DSURQE mouse atlas {cite}`Dorr2008-cc` - a high-resolution *ex vivo* structural MRI average of the mouse brain +accompanied with anatomical annotations provided by [the Mouse Imaging Centre (MICe)](https://www.mouseimaging.ca/) in Toronto, Canada. +The DSURQE template, mask and label files are all automatically pulled from the [original repository](https://www.mouseimaging.ca/repo/DSURQE_40micron/Dorr_2008_Steadman_2013_Ullmann_2013_Richards_2011_Qiu_2016_Egan_2015_40micron/) upon installing RABIES. + +## The EPICOMMON mouse template + +When using the `--bold_only` parameter to execute RABIES without structural scans, the common space +template file is changed for an alternative EPI reference template. This so-called EPICOMMON template +([available on zenodo](https://doi.org/10.5281/zenodo.19069284)) was generated by executing the initial +stages of the `--bold_only` RABIES preprocessing - 3D EPI generation, EPI inhomogeneity correction and unbiased +template generation from the EPI - using 5 independent fMRI datasets from different sites {cite}`Desrosiers-Gregoire2024-ou`. +The EPI template then inherited all mask and label files from the DSURQE atlas through non-linear registration. + +Providing an EPI template as target greatly improves the robustness of EPI registration, since image contrast and +shape are much more similar. + + +```{seealso} +[How to override the default common space template](../how_to/change_template.md) for +customising the input template and associated files. +``` diff --git a/docs/reference/workflows.md b/docs/reference/workflows.md new file mode 100644 index 00000000..b2b85fc3 --- /dev/null +++ b/docs/reference/workflows.md @@ -0,0 +1,119 @@ +# Workflow reference + +RABIES pipelines are built as [Nipype](https://nipype.readthedocs.io/en/latest/) +workflows: each node is a processing step, and the required inputs and outputs +define the links between nodes. + +The docstrings below are included directly from the RABIES source, so they +describe the version of the code these docs were built from. For what each +workflow is for and why it exists, see +[The preprocessing workflow](../explanation/preprocessing.md) and +[The confound correction workflow](../explanation/confound_correction.md). + +## Preprocessing + +(wf_bold_ref)= + +### `rabies.preprocess_pkg.bold_ref.init_bold_reference_wf` + +[Source](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/bold_ref.py) + +```{literalinclude} ../../rabies/preprocess_pkg/bold_ref.py +:start-after: gen_bold_ref_head_start +:end-before: gen_bold_ref_head_end +``` + +(wf_hmc)= + +### `rabies.preprocess_pkg.hmc.init_bold_hmc_wf` + +[Source](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/hmc.py) + +```{literalinclude} ../../rabies/preprocess_pkg/hmc.py +:start-after: hmc_wf_head_start +:end-before: hmc_wf_head_end +``` + +(wf_motion_params)= + +### `rabies.preprocess_pkg.hmc.EstimateMotionParams` + +[Source](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/hmc.py) + +```{literalinclude} ../../rabies/preprocess_pkg/hmc.py +:start-after: motion_param_head_start +:end-before: motion_param_head_end +``` + +(wf_inho_correction)= + +### `rabies.preprocess_pkg.inho_correction.init_inho_correction_wf` + +[Source](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/inho_correction.py) + +```{literalinclude} ../../rabies/preprocess_pkg/inho_correction.py +:start-after: inho_correction_head_start +:end-before: inho_correction_head_end +``` + +(wf_commonspace_reg)= + +### `rabies.preprocess_pkg.commonspace_reg.init_commonspace_reg_wf` + +[Source](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/commonspace_reg.py) + +```{literalinclude} ../../rabies/preprocess_pkg/commonspace_reg.py +:start-after: commonspace_wf_head_start +:end-before: commonspace_wf_head_end +``` + +(wf_cross_modal_reg)= + +### `rabies.preprocess_pkg.registration.init_cross_modal_reg_wf` + +[Source](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/registration.py) + +```{literalinclude} ../../rabies/preprocess_pkg/registration.py +:start-after: cross_modal_reg_head_start +:end-before: cross_modal_reg_head_end +``` + +(wf_bold_resampling)= + +### `rabies.preprocess_pkg.resampling.init_bold_preproc_trans_wf` + +[Source](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/resampling.py) + +```{literalinclude} ../../rabies/preprocess_pkg/resampling.py +:start-after: bold_resampling_head_start +:end-before: bold_resampling_head_end +``` + +(wf_mask_resampling)= + +### `rabies.preprocess_pkg.resampling.init_mask_preproc_trans_wf` + +[Source](https://github.com/CoBrALab/RABIES/blob/master/rabies/preprocess_pkg/resampling.py) + +```{literalinclude} ../../rabies/preprocess_pkg/resampling.py +:start-after: mask_resampling_head_start +:end-before: mask_resampling_head_end +``` + +## Confound correction + +(wf_confound_correction)= + +### `rabies.confound_correction_pkg.confound_correction.init_confound_correction_wf` + +[Source](https://github.com/CoBrALab/RABIES/blob/master/rabies/confound_correction_pkg/confound_correction.py) + +```{literalinclude} ../../rabies/confound_correction_pkg/confound_correction.py +:start-after: confound_wf_head_start +:end-before: confound_wf_head_end +``` + +```{seealso} +[How to contribute to RABIES](../how_to/contribute.md) covers writing a new +Nipype interface and connecting it into one of these workflows. +``` diff --git a/docs/requirements.txt b/docs/requirements.txt index 3e7f6c1e..4d546ec2 100644 --- a/docs/requirements.txt +++ b/docs/requirements.txt @@ -1,10 +1,17 @@ -sphinx==5.0 +# Docs toolchain. Upper bounds only where a real incompatibility exists, so +# routine upgrades do not need a change here. +# +# The Python version is capped at 3.12 in .readthedocs.yml because this file +# installs rabies (for sphinxcontrib-programoutput, which runs `rabies --help` +# at build time) and rabies requires Python <3.13. +sphinx myst-parser -sphinx-rtd-dark-mode sphinx-rtd-theme +sphinx-rtd-dark-mode +sphinx-design +sphinx-copybutton sphinxcontrib-bibtex sphinxcontrib-programoutput -jinja2==3.1.6 -pillow==12.3.0 +jinja2 +pillow rabies==0.6.1 -traits<7.0 diff --git a/docs/running_the_software.md b/docs/running_the_software.md deleted file mode 100644 index ba3a63c2..00000000 --- a/docs/running_the_software.md +++ /dev/null @@ -1,202 +0,0 @@ -# Running The Software - -## Input data in BIDS standard - -The input dataset must be organized according to the [BIDS data structure](https://bids.neuroimaging.io/){cite}`Gorgolewski2016-zm`. RABIES will iterate through all subjects found to contain a functional file (see section on BIDS filters below), and will also iterate according to sessions and runs found within each subject if available. If anatomical scans are used for preprocessing (i.e. not using `--bold_only`), each functional scan will be matched to one corresponding anatomical scan of the same subject/session (using BIDS filters for the anatomical image, see below). - -### Directory structure for an example dataset -* Our [example dataset](http://doi.org/10.5281/zenodo.8349029) has the following BIDS structure: - - - - - - - - - - - https://zenodo.org/record/8349029/preview/test_dataset.zip
- ├── sub-PHG001
- │ └── ses-3
- │     ├── anat
- │     │ ├── sub-PHG001_ses-3_acq-RARE_T2w.json
- │     │ └── sub-PHG001_ses-3_acq-RARE_T2w.nii.gz
- │     └── func
- │         ├── sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold.json
- │         └── sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold.nii.gz
- └── sub-PHG002
-     └── ses-3
-         ├── anat
-         │ ├── sub-PHG002_ses-3_acq-RARE_T2w.json
-         │ └── sub-PHG002_ses-3_acq-RARE_T2w.nii.gz
-         └── func
-             ├── sub-PHG002_ses-3_task-rest_acq-EPI_run-1_bold.json
-             └── sub-PHG002_ses-3_task-rest_acq-EPI_run-1_bold.nii.gz
-

-

-

- -8 directories, 8 files -

-

-
- - - -### BIDS filters to identify functional and structural images -By default, RABIES will use the 'bold' or 'cbv' suffix to identify functional scans and the 'T1w' or 'T2w' suffix for structural scans. Files which don't match the BIDS filters are ignored. However, the BIDS filters can also be customized with the `--bids_filter` parameter during the preprocessing stage. This can be useful for instance if the default is not enough to find the right set of scans. The custom BIDS filter must be formated into a JSON file with the functional filter under 'func' and structural filter under 'anat' (see example below for the default parameters): -```json -{ - "func": { - "suffix":["bold","cbv"] - }, - "anat": { - "suffix":["T1w","T2w"] - } -} -``` - -## Command Line Interface - -RABIES is executed using a command line interface, within a terminal. The software is divided into three main processing stages: preprocessing, confound correction and analysis. Accordingly, the command line interface allows for three different mode of execution, corresponding to the processing stages. So first, when executing the software, one of the processing stage must be selected. Below you can find the general --help message printed with `rabies --help`, which provides a summary of each processing stage together with options for parallel processing and memory management. Then, the --help associated to each processing stage, i.e. `preprocess`, `confound_correction` and `analysis`, describes in more detail the various parameters available to adapt image processing according to the user needs. Click on the corresponding --help to expand: - -
rabies --help -

- -```{program-output} rabies --help -``` - -

-
- -
rabies preprocess --help -

- -```{program-output} rabies preprocess --help -``` - -

-
- -
rabies confound_correction --help -

- -```{program-output} rabies confound_correction --help -``` - -

-
- - -
rabies analysis --help -

- -```{program-output} rabies analysis --help -``` - -

-
- - -## Example execution syntax -The following section provides examples describing the basic syntax for running the RABIES command line interface. - - -**preprocess** -```sh -rabies -p MultiProc preprocess input_BIDS/ preprocess_outputs/ --apply_STC --TR 1.2 --commonspace_reg masking=true,brain_extraction=false,template_registration=SyN,fast_commonspace=false -``` -First, we have to preprocess the dataset before it can be analyzed. In this example, we are running the RABIES preprocessing on the dataset found in the `input_BIDS/` folder, formatted according to the BIDS standard, and the outputs from RABIES are stored in the `preprocess_outputs/` folder. Additional execution parameters were specified: -* `-p MultiProc` will execute the pipeline in parallel using the local threads available. Notice that this parameter is specified before the processing stage, because it is one of the `Execution Options` affiliated to the `rabies --help`. -* `--apply_STC` is a boolean variable which, when selected, will apply slice timing correction during preprocessing, which is not applied by default in RABIES. -* `--TR 1.2` specifies the repetition time (TR) of the fMRI images that are processed, which must be defined to apply slice timing correction appropriately. Notice that this parameter must be provided with an argument, here `1.2` for TR = 1.2sec, and this is done by writing down the argument with a space dividing the associated parameter. -* `--commonspace_reg masking=true,brain_extraction=false,template_registration=SyN,fast_commonspace=false` this argument manages the options for the commonspace registration step. Some arguments, including `--commonspace_reg`, take multiple parameters as input, where each parameter-value pairs follow the syntax of `parameter=value`. In this case, we are using the masking optino with `masking=true`, which use available brain masks from inhomogeneity correction to drive the registration operations, and we specify a non-linear registration to the commonspace template with `template_registration=SyN`. - -**confound_correction** -```sh -rabies -p MultiProc confound_correction preprocess_outputs/ confound_correction_outputs/ --nuisance_regressors WM_signal CSF_signal vascular_signal mot_6 --smoothing_filter 0.3 -``` -Next, after completing preprocessing, in most cases the data should be corrected for potential confounds prior to analysis. This is done in the confound correction stage, where confounds are modelled and regressed from the data. In this example we correct the preprocessed data found in the `preprocess_outputs/` folder and store the cleaned outputs in the `confound_correction_outputs/` folder. Among the range of options available for confound correction, we define in this example three parameters: -* `--nuisance_regressors` is the option to regress nuisance timeseries from the data, i.e., confound regression. This parameter takes a list as input, where each argument in the list is seperated by a space as follow `WM_signal CSF_signal mot_6`. This list defines which nuisance timeseries are going to model confounds during confound regression, in this case, the WM and CSF mean signals together with the 6 rigid realignment parameters from head motion realignment. -* `--smoothing_filter` will additionally apply Gaussian spatial smoothing, where in this case, a filter size of `0.3` mm is specified. - -**analysis** -```sh -rabies -p MultiProc analysis confound_correction_outputs analysis_outputs/ --group_ica apply=true,dim=30,random_seed=1 -``` -Finally, after conducting preprocessing and confound correction, certain analyses can be run within RABIES. In this case, the cleaned outputs found in `confound_correction_outputs/` are going to be analyzed, with analysis outputs found in `analysis_outputs/`. We perform a group independent component analysis (ICA) with 30 components by providing `--group_ica apply=true,dim=30,random_seed=1` to the command. - -## Execution syntax with containerized installation (Apptainer and Docker) - -Containers are independent computing environments which have their own dependencies installed to ensure consistent and reliable -execution of the software across computing platforms. The main difference in execution syntax when running a container, as opposed to the examples above, is that the paths between the local environment where the data is stored must be 'linked' to the container's internal paths. All relevant directories containing data that will be used by RABIES must be related to a container internal path, and this is done using `-B` for Apptainer and `-v` for Docker. See below for examples: - -### Apptainer execution - -**preprocess** -```sh -apptainer run -B $PWD/input_BIDS:/input_BIDS:ro \ --B $PWD/preprocess_outputs:/preprocess_outputs/ \ -/path_to_apptainer_image/rabies.sif -p MultiProc preprocess /input_BIDS/ /preprocess_outputs/ --apply_STC --TR 1.2 --commonspace_reg masking=true,brain_extraction=false,template_registration=SyN,fast_commonspace=false -``` -Apptainer containers are stored in image files, for instance `rabies.sif`. `apptainer run /path_to_apptainer_image/rabies.sif` will execute the image, in this case the RABIES pipeline, and the same rules for the command line interface then apply as previously demonstrated. However, the container must gain access to the relevant folders for running RABIES, in this case an input folder and an output folder, and this is done with `-B`: -* `-B $PWD/input_BIDS:/input_BIDS:ro`: this argument relates the BIDS input folder found in `$PWD/input_BIDS` to an internal path to the container, which we call `/input_BIDS`. The inputs are thus accessed according to this path in the RABIES arguments with `/input_BIDS/`. the `:ro` means that the container is only provided reading permissions at this location. -* `-B $PWD/preprocess_outputs:/preprocess_outputs/`: same as with the `/input_BIDS/`, but now we are relating a desired output directory `$PWD/preprocess_outputs` to `/preprocess_outputs`, and the container has writing permissions at this path since `:ro` is not present. - - -**confound_correction** -```sh -apptainer run -B $PWD/input_BIDS:/input_BIDS:ro \ --B $PWD/preprocess_outputs:/preprocess_outputs/ \ --B $PWD/confound_correction_outputs:/confound_correction_outputs/ \ -/path_to_apptainer_image/rabies.sif -p MultiProc confound_correction /preprocess_outputs/ /confound_correction_outputs/ --nuisance_regressors WM_signal CSF_signal vascular_signal mot_6 --smoothing_filter 0.3 -``` -The required paths are similarly provided for the confound correction stage. Note here that the path to `$PWD/input_BIDS` is still linked to the container, even though it is not explicitely part of the arguments during the confound correction call. This is necessary since the paths used in the preprocessing steps still need to be accessed at later stages, and there will be an error if the paths are not kept consistent across processing steps. - -**analysis** -```sh -apptainer run -B $PWD/input_BIDS:/input_BIDS:ro \ --B $PWD/preprocess_outputs:/preprocess_outputs/ \ --B $PWD/confound_correction_outputs:/confound_correction_outputs/ \ --B $PWD/analysis_outputs:/analysis_outputs/ \ -/path_to_apptainer_image/rabies.sif -p MultiProc analysis /confound_correction_outputs /analysis_outputs/ --group_ica apply=true,dim=30,random_seed=1 -``` -The same logic applies at the analysis stage. -
- -### Docker execution -```sh -docker run -it --rm --user $(id -u) \ --v $PWD/input_BIDS:/input_BIDS:ro \ --v $PWD/preprocess_outputs:/preprocess_outputs/ \ -gabdesgreg/rabies:tagname -p MultiProc preprocess /input_BIDS/ /preprocess_outputs/ --apply_STC --TR 1.2 --commonspace_reg masking=true,brain_extraction=false,template_registration=SyN,fast_commonspace=false -``` -The syntax in Docker is very similar to Apptainer, except that `-B` is replaced by `-v`, and further parameters may be needed (e.g. `-it`, `--rm`). `--user $(id -u)` can be added to mitigate writing permission issues when using Docker. Note that 'tagname' should be replaced by the proper RABIES version you are using (e.g. 0.4.8). - - -## Additional Resources - -* **Workshop and tutorial for RABIES:** - * [Hands-on tutorial](https://github.com/CoBrALab/RABIES_tutorial) on RABIES (originally developed for and presented at the INCF Neuroinformatics Assembly 2023). - * A workshop providing a complete software overview was [recorded and posted online](https://www.youtube.com/watch?v=LZohKlUgycc&t=2766s&ab_channel=DouglasResearchCentre) on February 2023. -* Conversion from Bruker raw to Nifti formats can be handled with [BrkRaw](https://brkraw.github.io/) (consult [associated documentation](https://github.com/CoBrALab/documentation/wiki/bruker2nifti-conversion) from the CoBrALab) -* [CoBrALab recommendations](https://github.com/CoBrALab/documentation/wiki/Running-RABIES-on-niagara) for using compute canada. - diff --git a/docs/troubleshooting.md b/docs/troubleshooting.md deleted file mode 100644 index 14140181..00000000 --- a/docs/troubleshooting.md +++ /dev/null @@ -1,41 +0,0 @@ -# Troubleshooting - -This page provides guidance on common issues encountered when using RABIES. - -## Checking Image Orientation with ITK-SNAP - -RABIES assumes that input data are properly oriented according to the NIfTI standard (RAS+ orientation). -Incorrectly oriented images are a common source of processing failures and unexpected results. -Before reporting bugs or troubleshooting other issues, it is critical to verify that your images have the correct orientation. - -### Installing ITK-SNAP - -[ITK-SNAP](https://www.itksnap.org/pmwiki/pmwiki.php) is a free, open-source medical image viewer that properly displays NIfTI orientation information. -RABIES uses ANTs/ITK tools under the hood. -Incorrect image orientation is one of the most common causes of registration failures in RABIES. -Since RABIES expects RAS orientation (Right–Anterior–Superior), you should always verify your anatomical scans in ITK-SNAP before running the pipeline. - -### How to Check Orientation in ITK-SNAP - -1. **Open your image in ITK-SNAP** - - Go to **File → Open Main Image…** - - Load your NIfTI anatomical scan. - -2. **Verify anatomical orientation** - - ITK-SNAP shows three orthogonal views (Axial, Coronal, Sagittal). - - Ensure anatomical structures appear where you expect them (e.g., nose = anterior, top of head = superior). - -3. **Compare your scan to a correctly oriented atlas/template** - - Open your reference atlas (e.g., SIGMA, Fischer rat, etc.) in another ITK-SNAP window. - - Verify that structures appear in similar positions and that the orientation labels match. - -4. **Check the orientation labels** - - Inspect axes labels around each view (**R/L**, **A/P**, **S/I**). - - Confirm they correspond to the real anatomical directions in your scan. - - Move the cursor: the crosshair should move consistently across all views (e.g., dragging right corresponds to anatomical right). - -## See Also - -- For registration-specific troubleshooting, see [Registration Troubleshooting](nested_docs/registration_troubleshoot.md) -- For QC-related guidance, see [Preprocessing QC](preproc_QC.md) -- When reporting bugs, refer to the [issue template](https://github.com/CoBrALab/RABIES/blob/master/.github/ISSUE_TEMPLATE/standard-bug-report.md) diff --git a/docs/tutorials/first_run.md b/docs/tutorials/first_run.md new file mode 100644 index 00000000..ca512434 --- /dev/null +++ b/docs/tutorials/first_run.md @@ -0,0 +1,351 @@ +# Your first RABIES run + +In this tutorial you will take two raw mouse fMRI scans and carry them all the +way through the three RABIES stages: `preprocess`, `confound_correction` and `analysis`. +Along the way you will look at the quality control images RABIES produces, so +that by the end you will have seen every part of the software in action. + +You do not need to understand every option you type. The point is to build a +feel for the shape of a RABIES run. Everything you use here is explained +elsewhere, and this page links to the relevant pages as it goes. + +## What you will end up with + +Three output folders that mirror the three stages of the software: + +```text +rabies_tutorial/ +├── bids_inputs/ the data you download in step 2 +├── preprocess_outputs/ step 3 +├── confound_correction_outputs/ step 5 +└── analysis_outputs/ step 6 +``` + +```{admonition} How long this takes +:class: important + +Steps 1 and 2 take a few minutes. Step 3 is real image processing on real +brains, and is the long one — expect it to run for an hour or more, depending +on how many cores your machine has. Steps 5 and 6 take a couple of minutes +each. Start step 3 and come back to it. +``` + +## Step 1: Install the container + +RABIES has a lot of non-Python dependencies, so you will use the prebuilt +container. Install [Apptainer](https://apptainer.org/docs/user/main/quick_start.html), +then build the RABIES image: + +```sh +apptainer build rabies.sif docker://ghcr.io/cobralab/rabies:latest +``` + +This downloads roughly a gigabyte and takes a few minutes. When it finishes you +will have a single file, `rabies.sif`, which contains the entire software +environment. + +Check that it runs: + +```sh +apptainer run rabies.sif --help +``` + +You should see the RABIES usage message, listing the three processing stages: +`preprocess`, `confound_correction` and `analysis`. That is the structure of +the whole tutorial. + +```{tip} +If you are on macOS or Windows, or Apptainer is not available to you, use +Docker instead — see [How to install RABIES](../how_to/install.md). Every +`apptainer run ... -B src:dst` below becomes `docker run ... -v src:dst`. +``` + +## Step 2: Get the example data + +Make a working directory and download the two-subject example dataset: + +```sh +mkdir -p rabies_tutorial && cd rabies_tutorial +curl -L -o test_dataset.zip https://zenodo.org/records/8349029/files/test_dataset.zip +unzip test_dataset.zip -d bids_inputs +``` + +Look at what you got: + +```sh +find bids_inputs -name '*.nii.gz' +``` + +Two subjects, each with one anatomical scan and one resting-state functional +scan: + +```text +bids_inputs/ +├── sub-PHG001/ses-3/ +│ ├── anat/sub-PHG001_ses-3_acq-RARE_T2w.nii.gz +│ └── func/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold.nii.gz +└── sub-PHG002/ses-3/ + ├── anat/sub-PHG002_ses-3_acq-RARE_T2w.nii.gz + └── func/sub-PHG002_ses-3_task-rest_acq-EPI_run-1_bold.nii.gz +``` + +Notice the file names. The `sub-`, `ses-`, `task-` and `T2w`/`bold` parts are +not decoration — this is the [BIDS](https://bids.neuroimaging.io/) naming +standard, and it is how RABIES finds the scans and pairs each functional image +with the anatomical image from the same session. Your own data will need to +look like this too. See [Input data requirements](../reference/bids_inputs.md). + +## Step 3: Preprocess + +Create the output folder, then run the preprocessing stage: + +```sh +mkdir -p preprocess_outputs + +apptainer run \ + -B $PWD/bids_inputs:/inputs:ro \ + -B $PWD/preprocess_outputs:/outputs \ + rabies.sif -p MultiProc \ + preprocess /inputs /outputs \ + --anatomical_resampling 0.15x0.15x0.15 \ + --anat_inho_cor multiotsu=true +``` + +Breakdown of the command: + +- `-B $PWD/bids_inputs:/inputs:ro` makes your input folder visible inside the + container under the name `/inputs`, read-only. The container cannot see any + path you do not bind, which is why the RABIES arguments say `/inputs` and + `/outputs` rather than the paths on your machine. +- `-p MultiProc` runs the pipeline in parallel across the cores of your + machine. Leave it out and everything runs one step at a time. +- `--anatomical_resampling 0.15x0.15x0.15` commands to resample the images + to 0.15mm isotropic resolution for each registration operation. This + resolution is a decent heuristic in practice for a mouse brain to + derive 'good-enough' alignment for standard EPI images while speeding + things up for high resolution structural images. We do not recommend registration + in non-isotropic resolution. +- `--anat_inho_cor multiotsu=true` activates the multiotsu option for inhomogeneity + correction of anatomical images, which carries a more aggressive correction + (for most datasets this is not necessary, you should not use this by default, + here it's because this particular dataset has an especially sharp intensity gradient). + +You did not specify an atlas. The container comes with the DSURQE mouse atlas +already installed, and RABIES uses it by default (see [built-in template files](../reference/template_files.md)). + +Let it run. You will see a stream of Nipype log lines as each node in the +workflow completes. + +## Step 4: Look at the registration quality control report + +This is a **crucial step** in practice — this is where you'll catch alignment issues that can +completely break later stages. Open the report folder: + +```sh +ls preprocess_outputs/preprocess_QC_report/ +``` + +Each subfolder holds PNG images for one key preprocessing step that requires visual +quality control (formally documented [here](../reference/qc_outputs.md)). +You will go through each subfolder in order, as later steps depend on the success +of the first ones. + +1. First open `preprocess_outputs/preprocess_QC_report/anat_inho_cor` which displays +the inhomogeneity correction step for the structural scans. + +```{figure} ../pics/sub-PHG001_ses-3_acq-RARE_T2w_inho_cor.png +:alt: Four-column figure showing the stages of structural inhomogeneity correction + +Four panels of the inhomogeneity correction report, representing the initial and +final iterations of inhomogeneity before/after brain masking. +``` +Your outputs should resemble the above example, where brain masking was successfully +delineating the brain edges, and the prominent intensity gradient in the image is +properly corrected for (partly thanks to the `multiotsu=true` option we selected +above). + +The success of this inhomogeneity correction is crucial for the registration step +below, as errors here can leave intensity biases that do not match across images, +or even cropped out brains from masking errors. + + +2. Then, you can consult the alignment of each image to the dataset average (of only +2 scans in this case) within `preprocess_QC_report/commonspace_reg_wf.Anat2Unbiased/`. + +```{figure} ../pics/sub-PHG001_ses-3_acq-RARE_T2w_RAS_inho_cor_registration.png +:alt: Overlap to the study template + +The alignment of your data (top) with the study template (bottom). The +outlines should follow the same anatomy. +``` +The alignment here determines whether an entire scanning session is properly aligned +to common space. An error here will result in bad anatomical masks that impact EPI +preprocessing steps below, and ultimately misaligned timeseries at the end of preprocessing. + +3. Next, you should validate the registration of the study template to the +external commonspace template in `preprocess_QC_report/commonspace_reg_wf.Unbiased2Atlas/`. + +```{figure} ../pics/test_dataset_atlas_registration.png +:alt: Overlap of the study template on the reference atlas + +The alignment of your study template (top) with the external commonspace template (bottom). The +outlines should follow the same anatomy. +``` +The alignment to the external atlas is the last step that can link up the set of brain masks that are +then used for inhomogeneity correction/registration of the EPI below, and also regulates the +eventual resampling of timeseries to common space. + +4. We then move on to the EPI preprocessing steps. You can open the `preprocess_QC_report/bold_inho_cor/` +folder, where you will validate the inhomogeneity correction of each EPI scan. + +```{figure} ../pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_inho_cor.png +:alt: Four-column figure showing the stages of functional inhomogeneity correction + +Same four panels of the inhomogeneity correction report as with the structural scan, +but now for the EPI scan. +``` +Again here you should find adequate brain masking and intensity correction. An error here +would impact the EPI registration below. + +5. Finally, open `preprocess_QC_report/EPI2Anat/`, which shows each functional scan +aligned onto its own anatomical scan: + +```{figure} ../pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_registration.png +:alt: Overlap of a functional scan on its anatomical scan + +The functional image (top) matched to the anatomical image (bottom) from the +same session. This registration is also what corrects susceptibility +distortion. +``` +It is unlikely that this registration step will work if **any** of the previous +steps failed significantly for a given run/session. + +On this example dataset, each preprocessing step should run without error. +When they do not on your own data, [How to troubleshoot registration](../how_to/troubleshoot_registration.md) +will help you fix it. + +## Step 5: Correct confounds + +Preprocessed data is not yet ready to analyse: the next stage is critical for +removing head motion and physiological artefacts that produce correlations +that corrupt functional connectivity. + +```sh +mkdir -p confound_correction_outputs + +apptainer run \ + -B $PWD/bids_inputs:/inputs:ro \ + -B $PWD/preprocess_outputs:/outputs \ + -B $PWD/confound_correction_outputs:/cc_outputs \ + rabies.sif -p MultiProc \ + confound_correction /outputs /cc_outputs \ + --frame_censoring FD_censoring=true,FD_threshold=0.05 \ + --nuisance_regressors mot_6 aCompCor_5 \ + --smoothing_filter 0.3 +``` + +Note that you bind `bids_inputs` again even though it is not in the RABIES +arguments. Each stage reads the file paths recorded by the previous one, so +every path used in step 3 has to stay reachable, at the same location, for the +rest of the pipeline. + +The options you passed: + +- `--frame_censoring FD_censoring=true,FD_threshold=0.05` applies censoring + of high-motion frames using a 0.05 mm threshold on framewise displacement. +- `--nuisance_regressors mot_6 aCompCor_5` models the signal using + the six head motion parameters plus the first 5 aCompCor principal + components derived from the combined WM-CSF masks, and subtracts + the variance explained by those regressors. +- `--smoothing_filter 0.3` applies 0.3 mm Gaussian spatial smoothing. + +That is a relatively modest correction. Choosing a correction strategy for a +real dataset is its own task, covered in +[How to optimise your confound correction strategy](../how_to/optimise_confound_correction.md). + +The cleaned timeseries land in +`confound_correction_outputs/confound_correction_datasink/cleaned_timeseries/`. + +## Step 6: Analyse + +Now compute seed-based connectivity using a pre-built seed, and generate data +quality reports using `--data_diagnosis`: + +```sh +mkdir -p analysis_outputs + +apptainer run \ + -B $PWD/bids_inputs:/inputs:ro \ + -B $PWD/preprocess_outputs:/outputs \ + -B $PWD/confound_correction_outputs:/cc_outputs \ + -B $PWD/analysis_outputs:/analysis_outputs \ + rabies.sif -p MultiProc \ + analysis /cc_outputs /analysis_outputs \ + --data_diagnosis \ + --seed_list SS_frontal_seed +``` + +Breakdown of the command: +- `--seed_list SS_frontal_seed` will use the pre-built seed in the frontal cortex + for deriving seed-based connectivity in each scan. +- `--data_diagnosis` will generate a set of data quality reports associated to the + functional connectivity analysis results. + +You have now carried out a basic functional connectivity analysis, and +generated the [scan-level spatiotemporal diagnosis report](../explanation/scan_diagnosis.md) +for each scan, which you can now open from `analysis_outputs/data_diagnosis_datasink/figure_temporal_diagnosis/` +and `analysis_outputs/data_diagnosis_datasink/figure_spatial_diagnosis/` subfolders. +This report is part of a larger [data quality assessment](../explanation/data_quality.md) +framework developed for rodent fMRI along with RABIES. Some group-level reports +were not generated here, as those require at minimum 3 scans. +It takes time to learn how to read these reports, and some experience to +interpret them accurately. For the sake of the tutorial, we will simply +validate whether we can observe a brain network from seed-based connectivity. Let's open the spatial +diagnosis report under `analysis_outputs/data_diagnosis_datasink/figure_spatial_diagnosis/`: + +```{figure} ../pics/sub-PHG001_ses-3_task-rest_acq-EPI_run-1_bold_spatial_diagnosis.png +:alt: Tutorial spatial diagnosis report + +Spatial diagnosis report for the first subject. +``` +In this report, the resulting seed-based connectivity can be visualised in the +last row labelled 'SBC network 0'. We can see in those two example scans that +the seed revealed a bilateral correlation structure that corresponds to +the somatomotor network anatomy of the mouse brain. This confirms that +this network was adequately mapped and measured (although more subtle +confounding effects could still exist, but that is a question beyond +the scope of the tutorial). + +These connectivity maps are saved and can be accessed as NIfTI files +in the `analysis_outputs/commonspace_analysis_datasink/seed_correlation_maps` +output folder, and could be fed into downstream statistical analyses in an +actual experiment. + +## What you have done + +You have run all three RABIES stages and produced a connectivity estimate from +raw scans. Specifically, you: + +- built the container and confirmed it runs +- preprocessed two subjects, registering them into a common atlas space +- inspected the registration quality control images +- removed motion and physiological confounds from the timeseries +- computed a seed-based connectivity on each scan +- confirmed the mapping of the somatomotor network using the scan-level + diagnostic report + +## Where to go next + +For your own data, the decisions the tutorial made for you become yours to +make: + +- Your data has to be in BIDS format — see + [Input data requirements](../reference/bids_inputs.md). +- Your registrations will not always succeed on the first try — see + [How to troubleshoot registration](../how_to/troubleshoot_registration.md). +- Confound correction should be tuned to your data, not copied from a tutorial + — see [How to optimise your confound correction strategy](../how_to/optimise_confound_correction.md) + and [How to assess data quality](../how_to/assess_data_quality.md). +- To understand what the preprocessing actually did, read + [The preprocessing workflow](../explanation/preprocessing.md). + diff --git a/docs/tutorials/index.md b/docs/tutorials/index.md new file mode 100644 index 00000000..a1f9501f --- /dev/null +++ b/docs/tutorials/index.md @@ -0,0 +1,23 @@ +# Tutorial + +A tutorial is a lesson. It takes you through a complete piece of work from +beginning to end so that you learn what RABIES does by doing it, on data that +is known to work. You are not expected to make any decisions along the way, and +you do not need to understand every step as it happens. + +```{toctree} +--- +maxdepth: 1 +--- +first_run +``` + +Once you have finished, the [how-to guides](../how_to/index.md) will help you +apply RABIES to your own data. + +```{seealso} +The CoBrALab also maintains a longer +[hands-on workshop](https://github.com/CoBrALab/RABIES_tutorial), originally +developed for the INCF Neuroinformatics Assembly 2023, and a +[recorded software overview](https://www.youtube.com/watch?v=LZohKlUgycc&t=2766s&ab_channel=DouglasResearchCentre). +``` diff --git a/rabies/parser.py b/rabies/parser.py index fe6526b8..c15b1d8a 100644 --- a/rabies/parser.py +++ b/rabies/parser.py @@ -580,7 +580,7 @@ def get_parser(): description= "Specify commonspace template and associated mask files. By default, RABIES\n" "provides the mouse DSURQE atlas\n" - "https://wiki.mouseimaging.ca/display/MICePub/Mouse+Brain+Atlases.\n" + "https://www.mouseimaging.ca/repo/DSURQE_40micron/Dorr_2008_Steadman_2013_Ullmann_2013_Richards_2011_Qiu_2016_Egan_2015_40micron/.\n" ) g_commonspace.add_argument( '--anat_template', action='store', type=Path,