--- jupytext: formats: ipynb,md:myst text_representation: extension: .md format_name: myst format_version: 0.13 jupytext_version: 1.11.5 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- ```{code-cell} ipython3 :tags: [remove-cell] from ndslib.config import jupyter_startup jupyter_startup() ``` (nibabel)= # Reading neuroimaging data with `Nibabel` The first step in almost any neuroimaging data analysis is reading neuroimaging data from files. We will start dipping our toes into how to do this with the nibabel software library. Many different file formats store neuroimaging data, but the NIfTI format (The acronym stands for Neuroimaging Informatics Technology Initiative, which was the group at the US National Institutes of Health that originally designed the file format) has gained significant traction in recent years. For example, it is the format required by the BIDS standard for storing raw MRI data. Because this file format is so common, we will focus on this format here and we will not go into detail on the other file formats that are supported by nibabel. Getting acquainted with nibabel and its operations will also allow us to work through some issues that arise when analyzing data from neuroimaging experiments. In the next section ({numref}`nibabel_transforms`), we will dive even deeper into the challenge of aligning different measurements to each other. Let's start by downloading some neuroimaging files and seeing how nibabel works. We'll use the data downloading functionality implemented in ndslib to get data from the [OpenNeuro archive](https://openneuro.org/). These two files are part of a study on the effects of learning on brain representations {cite}`Beukema2019-ft`. From this dataset, we will download just one run of a functional MRI experiment and a corresponding T1-weighted image acquired for the same subject in the same session. ```{code-cell} ipython3 from ndslib.data import download_bids_dataset download_bids_dataset() ``` After running this code, we should have a valid BIDS dataset inside a folder called `ds001233`. If you examine your filesystem, you will see that the files have been named and organized according to the BIDS standard. This will come in handy in a little bit. ``` ds001233 ├── dataset_description.json └── sub-17 └── ses-pre ├── anat │   ├── sub-17_ses-pre_T1w.json │   └── sub-17_ses-pre_T1w.nii.gz └── func ├── sub-17_ses-pre_task-cuedSFM_run-01_bold.json └── sub-17_ses-pre_task-cuedSFM_run-01_bold.nii.gz ``` In the meanwhile, let's start using nibabel to access the data: ```{eval-rst} .. index:: single: Nibabel ``` ```{code-cell} ipython3 import nibabel as nib img_bold = nib.load("ds001233/sub-17/ses-pre/func/sub-17_ses-pre_task-cuedSFM_run-01_bold.nii.gz") ``` What kind of object is `img_bold`? ```{code-cell} ipython3 type(img_bold) ``` As you can see, it's a `Nifti1Image` object. This is a special nibabel object for handling NIfTI files. One of the principles that guide the design of nibabel is that operations that take time, or that require a lot of the computer's memory are delayed and only executed when there is no choice. Instead, nibabel starts by reading a relatively small section of the file, called the file's "header" that contains only meta-data: information about the data, rather than the data itself. For example, the header contains information about the size of the data array that is stored inside the file. So, when we called `nibabel.load` above none of the actual data was read. Instead, nibabel reads the header information about this data and could be prepared to also give you the data if and when you want it. This is useful for cases where we need to know some property of the data, such as the size of the data array, but we don't yet need to read any of the data itself. For example: ```{code-cell} ipython3 print(img_bold.shape) ``` This tells us that the data in this file is a four-dimensional array with dimensions of 96 by 96 by 66 voxels, collected in 241 time points. There is other information stored in the header, and we'll look at some of this information a little bit later. First, we will use the object's `get_fdata` method to read the data from the file into memory. The 'f' in this method name indicates to us that nibabel is going to take the data that is stored in the file and do its best to represent it as floating point numbers, regardless of the way that the data is stored in the file. For example ```{code-cell} ipython3 data_bold = img_bold.get_fdata() type(data_bold) ``` Now it's a regular Numpy array. Everything that applies to the arrays that you already saw in previous chapters applies to this array. You can compute on it, visualize it, and so forth. What data type does it use? ```{code-cell} ipython3 print(data_bold.dtype) ``` As promised, the data, regardless of how it was stored in the file, is given to us in the 64-bit floating point Numpy data type. There are other ways to read data from nibabel, that allow you more fine-grained control of how much of the data you are reading, how and whether you store it in memory, and so on, but we will not go into that here. We recommend reading the nibabel documentation to get acquainted with these details. One of the reasons that we don't think this is worth expounding upon is that the `get_fdata` does the right thing for you in 99% of the cases that we have encountered ourselves. So, this provides a simple way to read the data and know what you are getting. What shape does this array have? ```{code-cell} ipython3 print(data_bold.shape) ``` It's exactly the shape you would expect from the information provided in the `Nifti1Image` that we saw above. There are many things we can do with this data, but one of the first things that you might want to do with a new dataset is to assess the quality of the data. This will help guide our analysis strategy and allow us to identify samples (time points, runs, subjects) that should be handled carefully (for example, excluded from further analysis, or considered with a grain of salt). +++ ### Assessing MRI data quality MRI data quality assessments are important because MRI data can vary quite a bit in terms of its quality and the conclusions from data analysis may be misguided if data quality is not taken into account. For example, the quality of measurements from young children and older adults could be affected by the tendency of participants from these populations to move more in the MRI scanner. There are many different ways to assess MRI data quality. We will start with a commonly-used metric that is also relatively easy to understand, by computing the temporal signal-to-noise ratio (or tSNR) of the measurement in each voxel. This is quantified as the mean signal in the voxel across time points, divided by the standard deviation of the signal across time points. To demonstrate the various concepts that we discussed above, we will use three different ways to compute this value. The first uses Numpy directly on the data that we read from the file: ```{eval-rst} .. index:: single: Temporal SNR ``` ```{code-cell} ipython3 import numpy as np tsnr_numpy = np.mean(data_bold, -1) / np.std(data_bold, -1) ``` The second way to compute relies on the Nipype software. We will not go into much depth about how to use the Nipype library, recommending instead that you refer to one of the resources provided below, but briefly: Nipype interfaces are objects that provide direct Python control of a variety of software, including functionality that implemented in many non-Python libraries, such as AFNI and FSL. Here, we are going to use a Nipype interface for a relatively straightforward Python-based calculation of SNR, very similar to the calculation we just did using Numpy. The way this is done here is by defining the interface object `TSNR()`, associating this object with an input file, the file that contains the BOLD data, and then calling the `.run()` method. This creates a file -- `tsnr.nii.gz` -- which contains the results of the computation. ```{eval-rst} .. index:: single: Nipype ``` ```{code-cell} ipython3 from nipype.algorithms.confounds import TSNR tsnr = TSNR() tsnr.inputs.in_file = 'ds001233/sub-17/ses-pre/func/sub-17_ses-pre_task-cuedSFM_run-01_bold.nii.gz' res = tsnr.run() ``` The output file ("tsnr.nii.gz") can be immediately loaded with nibabel: ```{code-cell} ipython3 tsnr_img = nib.load("tsnr.nii.gz") tsnr_data = tsnr_img.get_fdata() ``` (nibabel_ex1)= #### Exercise Are the results from the Numpy calculation and the Nipype calculation identical? What happens if you exclude all the voxels for which the standard deviation is zero or very small (e.g., smaller than 0.001)? +++ ```{eval-rst} .. index:: single: BIDS apps ``` ```{eval-rst} .. index:: single: MRIQC ``` A third option to calculate quality assessments for MRI data will use a BIDS app. The MRIQC app, already mentioned before in {numref}`bids` is implemented as a Docker container. It takes a valid BIDS dataset as input and produces a range of quality assessments as output. Let's see how this works in practice. In a terminal prompt, we first create a folder for the outputs of processing: ``` $ mkdir output ``` Then, we use Docker -- much as we did in section {numref}`docker` -- to point the MRIQC Docker container to the BIDS dataset: ``` $ docker run -it --rm -v $(pwd)/ds001233:/data:ro -v $(pwd)/output:/out nipreps/mriqc:latest /data /out participant --participant_label 17 ``` Because the dataset is rather small: one fMRI run and one T1-weighted scan, it takes just a few minutes to run. When this is completed, we should have the following contents in the 'output' folder: ``` output ├── dataset_description.json ├── logs ├── sub-17 │   └── ses-pre │   ├── anat │   │   └── sub-17_ses-pre_T1w.json │   └── func │   └── sub-17_ses-pre_task-cuedSFM_run-01_bold.json ├── sub-17_ses-pre_T1w.html └── sub-17_ses-pre_task-cuedSFM_run-01_bold.html ``` Notice that the outputs of this processing are also organized as a (derivatives) BIDS dataset. There are no neuroimaging data files in here, but the files named `sub-17_ses-pre_T1w.json` and `sub-17_ses-pre_task-cuedSFM_run-01_bold.json` each contain a dictionary of quality assessment metrics, including for `sub-17_ses-pre_task-cuedSFM_run-01_bold.json`, a field with the median tSNR value that was calculated after applying motion correction, in addition to many other metrics. In addition, the html files saved at the top level of the outputs contain some visual diagnostics of image quality. This example again demonstrates the importance of the BIDS standard. Without it, the dockerized MRIQC BIDS app would not know how to find and understand the data that is stored in the input folder, and we wouldn't necessarily have an immediate intuitive understanding of the outputs that were provided. (nibabel_ex2)= #### Exercise Is the MRIQC median tSNR larger or smaller than the ones calculated with Nipype and Numpy? Why is that? ## Additional resources There is extensive literature on approaches to quality assessment of MRI data. The [MRIQC documentation](https://mriqc.readthedocs.io/) is a good place to start reading about the variety of metrics that can be used to automatically characterize image quality. In particular, to see some of the impacts of MRI quality assessments on the conclusions from MRI studies, you can look at some of the original research literature on this topic. Across different MRI modalities, researchers have shown that conclusions from MRI studies can be confounded particularly by the impact that motion has on MRI data quality. This has been identified in diffusion MRI {cite}`Yendiki2014-be`, resting state functional connectivity {cite}`Power2012-nf` and in the analysis of structural MRI {cite}`Reuter2015-jz`.