Commit 8ffbcef9 authored by William Clarke's avatar William Clarke
Browse files

Update documentation for the new basis tools for MM. Touch up notebooks and...

Update documentation for the new basis tools for MM. Touch up notebooks and add one for conversion of LCModel basline. Check uncertianty validation.
parent ee1941b7
Pipeline #14859 canceled with stages
......@@ -42,4 +42,13 @@ scale
diff
****
| *Example* :code:`basis_tools diff --add_or_sub sub mega_on mega_off mega_diff`
| Form a basis set for a difference method using two other basis set. Add or subtract using :code:`--add_or_sub {'add'|'sub}`.
\ No newline at end of file
| Form a basis set for a difference method using two other basis set. Add or subtract using :code:`--add_or_sub {'add'|'sub}`.
add_set
*******
| *Example* :code:`basis_tools add_set --add_MM basis_without_mm/ basis_with_default_mm/`
| Add a (predefined) set of Gaussian peaks to a basis set. Three defualt sets are defined:
| 1) The FSL 'default' with peaks at 0.9, 1.2, 1.4, 1.7 ppm and a linked set at 2.08 & 3.0 ppm
| 2) The (experimental) MEGA edited MM peaks at 0.915 & 3.0 ppm (ratio of 3.75:2.0).
| 3) A water peak at 4.65 ppm
| The other options (:code:`--gamma --sigma`) allow a custom set of peaks to be specified.
......@@ -162,10 +162,6 @@ Below are detailed explanations of some of the optional arguments in the wrapper
Polynomial baseline order. Set to -1 to remove the baseline altogether.
:code:`--metab_groups`
Group metabolites into sub-groups that get their own lineshape parameters (shift and broadening). This can either be a list of integers (one per metabolite) from 0 to the max number of groups minus one. Or it could be a list of metabolites to be grouped. E.g. using the flag :code:`--metab_groups Mac NAA+NAAG+Cr` then the Mac spectrum will have its own group, the NAA, NAAG, and Cr will be in a different group, and all other metabolites in a 3rd group. Other possibilities are combine_all and separate_all, where metabs are combined into a single group or separated into distinct groups respectively.
:code:`--add_MM`
Add macromolecule peaks at the following frequencies: 0.9, 1.2, 1.4, 1.7 ppm and a doublet at 2.08 & 3.0 ppm
:code:`--add_MM_MEGA`
Add linked macromolecule peaks at 0.915 & 3.0 ppm (ratio of 3.75:2.0). This option is experimental!
:code:`--lorentzian`
By default the lineshape is a Voigt (lorentizian+gaussian). Use this flag to set to Lorentzian.
:code:`--ind_scale`
......@@ -189,7 +185,6 @@ The wrapper scripts can also take a configuration file as an input. For example,
ppmlim = [0.3,4.1]
metab_groups = combine_all
TE = 11
add_MM
report
The the following calls to :code:`fsl_mrs` or :code:`fsl_mrsi` are equivalent:
......@@ -199,7 +194,7 @@ The the following calls to :code:`fsl_mrs` or :code:`fsl_mrsi` are equivalent:
::
fsl_mrs --ppmlim .3 4.1 --metab_groups combine_all --TE 11 --add_MM --report
fsl_mrs --ppmlim .3 4.1 --metab_groups combine_all --TE 11 --report
......
......@@ -13,7 +13,7 @@ For an in depth discussion of the effects of MM basis spectra choice on fitting
Synthetic MM
~~~~~~~~~~~~
Syntheic MM basis spectra may be added using the :code:`--add_MM` flag with :code:`fsl_mrs` or :code:`fsl_mrsi`. In the interactive environment the same can be achieved by calling the method :code:`add_MM_peaks` of :code:`fsl_mrs.core.MRS`.
Synthetic MM basis spectra can be added to a basis set using :code:`basis_tools add_set --add_MM basis_in basis_out`. In the interactive environment the same can be achieved by calling the method :code:`add_default_MM_peaks` in a basis object (of type :code:`fsl_mrs.core.Basis`).
By default this option will add the following basis spectra (in separate metabolite groups) to the basis sets.
......@@ -27,7 +27,9 @@ By default this option will add the following basis spectra (in separate metabol
MM17, 1.7, 2, 10/20
MM21, "2.08, 2.25, 1.95, 3.0", "1.33, 0.22, 0.33, 0.4", 10/20
Additional peaks may be added int he interactive environment by calling :code:`add_MM_peaks` with optional arguments to override the defaults.
Additional peaks may be added in using :code:`basis_tools add_set --gamma ...`, or in the interactive environment by calling :code:`basis.add_MM_peaks`.
Note that when fitting using the synthetic macromolecular peaks the :code:`--metab_groups` option should be used to assign each synthetic MM peak to its own group. E.g., :code:`--metab_groups MM09 MM12 MM14 MM17 MM21`
References
~~~~~~~~~~
......
......@@ -61,9 +61,9 @@
"source": [
"## 3. Fitting\n",
"\n",
"The below is a call to the main MRSI wrapper script. Note the use of the `--add_MM` flag as this metabolite basis does not contain a macromolecule basis. Also note that the mrsi dataset provided is accompagnied by a JSON sidecar file which contains some useful information such as the echo time, the central frequency, and the dwell time.\n",
"The below is a call to the main MRSI wrapper script. Note that the basis set used here contains FSL-MRS default MM, for these to work properly we must assign them to individual metabolite groups. This is done by providing the `--metab_groups MM09 MM12 MM14 MM17 MM21` argument. To see how we created the basis set have a look at `example_usage/Example basis spectra conversion.ipynb`. \n",
"\n",
"The script will by default run in parallel on the available CPU cores. Depending on hardware this should take a few minutes.\n"
"The script will by default run in parallel on the available CPU cores. Depending on hardware this should take a few minutes."
]
},
{
......@@ -73,16 +73,17 @@
"outputs": [],
"source": [
"%sx fsl_mrsi --data example_data/example_mrsi/mrsi.nii.gz \\\n",
" --basis example_data/example_mrsi/3T_slaser_32vespa_1250.BASIS \\\n",
" --basis example_data/example_mrsi/3T_slaser_32vespa_1250_noTMS_defaultMM \\\n",
" --output MRSI/example_mrsi_fit \\\n",
" --mask example_data/example_mrsi/mask.nii.gz \\\n",
" --h2o example_data/example_mrsi/wref.nii.gz \\\n",
" --tissue_frac MRSI/mrsi_seg_wm.nii.gz MRSI/mrsi_seg_gm.nii.gz MRSI/mrsi_seg_csf.nii.gz \\\n",
" --add_MM \\\n",
" --baseline_order 2 \\\n",
" --metab_groups MM09 MM12 MM14 MM17 MM21\\\n",
" --combine PCho GPC --combine Cr PCr --combine NAA NAAG --combine Glu Gln --combine Glc Tau \\\n",
" --ignore Gly Gly_1 HG \\\n",
" --overwrite"
" --ignore Gly \\\n",
" --overwrite \\\n",
" --report"
]
},
{
......@@ -135,7 +136,7 @@
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"display_name": "Python 3.8.12 ('fsl_mrs')",
"language": "python",
"name": "python3"
},
......@@ -149,7 +150,12 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.10"
"version": "3.8.12"
},
"vscode": {
"interpreter": {
"hash": "6c11d37c82810953c5a08a185ec224dab920e965fab2a4fd7bf60d843c04e747"
}
}
},
"nbformat": 4,
......
%% Cell type:markdown id: tags:
# Example of MRSI fitting on the command line
This notebook demos the process of fitting an MRSI scan using the command line scripts included in FSL-MRS.
### Contents:
- [1. Reconstruction, Processing and Data Conversion](#1.-Reconstruction,-processing-and-data-conversion)
- [2. Segmentation](#2.-Tissue-segmentation)
- [3. Fitting](#3.-Fitting)
- [4. Visualisation of fit](#4.-Visualisation-of-fit)
Will Clarke
June 2020
University of Oxford
%% Cell type:markdown id: tags:
## 1. Reconstruction, processing and data conversion
MRSI reconstruction (from k-space data) can be highly specialised depending on the sequence. For example reconstruction of non-cartesian trajectories (e.g. spirals or concentric-rings) requires regridding or use of the NUFFT. Due to the specialised nature of MRSI reconstruction FSL-MRS does not contain tools for this step.
Simmilarly post-processing of MRSI data is commonly done as part of the reconstruction process. Though many of the processing tools in FSL-MRS can be run on MRSI data they have not been created with MRSI in mind.
This example therefore assumes that you are able to provide reconstructed and processed data ready for fitting. Data can be converted to NIfTI from the Siemens DICOM format using spec2nii. The authors of FSL-MRS and spec2nii are happy to add supported formats if example data and interpretation is provided.
%% Cell type:markdown id: tags:
## 2. Tissue segmentation
Tissue segmentation is required to have meaningful scaling of water referenced metabolite concentrations. mrsi_segment will produce three files, corresponding to the GM, WM and CSF FAST tissue segmentations registered to the MRSI volume.
Run tissue segmentation on the packaged T1 data and mask using the MRSI data file. Here we provide a (partial) .anat file produced by [fsl_anat](https://fsl.fmrib.ox.ac.uk/fsl/fslwiki/fsl_anat) to speed up execution.
This step requires an installation of FSL compatible with fslpy.
%% Cell type:code id: tags:
``` python
#%%capture
%sx mkdir -p MRSI
%sx mrsi_segment -o MRSI -a example_data/example_mrsi/T1.anat example_data/example_mrsi/mrsi.nii.gz
```
%% Cell type:markdown id: tags:
## 3. Fitting
The below is a call to the main MRSI wrapper script. Note the use of the `--add_MM` flag as this metabolite basis does not contain a macromolecule basis. Also note that the mrsi dataset provided is accompagnied by a JSON sidecar file which contains some useful information such as the echo time, the central frequency, and the dwell time.
The below is a call to the main MRSI wrapper script. Note that the basis set used here contains FSL-MRS default MM, for these to work properly we must assign them to individual metabolite groups. This is done by providing the `--metab_groups MM09 MM12 MM14 MM17 MM21` argument. To see how we created the basis set have a look at `example_usage/Example basis spectra conversion.ipynb`.
The script will by default run in parallel on the available CPU cores. Depending on hardware this should take a few minutes.
%% Cell type:code id: tags:
``` python
%sx fsl_mrsi --data example_data/example_mrsi/mrsi.nii.gz \
--basis example_data/example_mrsi/3T_slaser_32vespa_1250.BASIS \
--basis example_data/example_mrsi/3T_slaser_32vespa_1250_noTMS_defaultMM \
--output MRSI/example_mrsi_fit \
--mask example_data/example_mrsi/mask.nii.gz \
--h2o example_data/example_mrsi/wref.nii.gz \
--tissue_frac MRSI/mrsi_seg_wm.nii.gz MRSI/mrsi_seg_gm.nii.gz MRSI/mrsi_seg_csf.nii.gz \
--add_MM \
--baseline_order 2 \
--metab_groups MM09 MM12 MM14 MM17 MM21\
--combine PCho GPC --combine Cr PCr --combine NAA NAAG --combine Glu Gln --combine Glc Tau \
--ignore Gly Gly_1 HG \
--overwrite
--ignore Gly \
--overwrite \
--report
```
%% Cell type:markdown id: tags:
## 4. Visualisation of fit
Now take a look at the outputs. A PNG of the fit to the average of all voxels is provided for a quick sanity check. The folders contain the following:
- concs : concentration for each metabolite or combined metabolites (subfolders contain different types of referencing)
- fit : the model prediction FID, the residual FID, and the baseline (also in the time domain).
- qc : QC parameters split per metabolite
- uncertainties : the std of the fit per metabolite
%% Cell type:code id: tags:
``` python
%sx ls -l MRSI/example_mrsi_fit
```
%% Cell type:markdown id: tags:
Now run the command below in a terminal in order to load the data in FSLeyes. You will need to install the NIfTI-MRS plugin for FSLeyes. Instructions for installation are available [online](https://open.win.ox.ac.uk/pages/wclarke/fsleyes-plugin-mrs/install.html).
Then follow the the instructions [here](https://open.win.ox.ac.uk/pages/wclarke/fsleyes-plugin-mrs/mrsi_results.html) to set it up in such a way that you can explore each voxel's individual fit.
To get started, after installing the plugin, run the following command:
```
fsleyes -smrs example_data/example_mrsi/T1.anat/T1.nii.gz example_data/example_mrsi/mrsi.nii.gz
```
%% Cell type:code id: tags:
``` python
%sx fsleyes -smrs example_data/example_mrsi/T1.anat/T1.nii.gz example_data/example_mrsi/mrsi.nii.gz
```
......
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Command-line conversion of LCModel basis set\n",
"This notebook demonstrates the conversion of an existing LCModel basis set (`.basis` format).\n",
"\n",
"In addition to the basic conversion step, three extra steps are performed:\n",
"1. Removal of the TMS peak from all spectra\n",
"2. Removal of duplicated and unwanted basis spectra\n",
"3. Addition of FSL-MRS default MM peaks\n",
"\n",
"This converted spectrum will subsequently be used to fit the example MRSI data.\n",
"\n",
"## Step 1 - conversion\n",
"We use the `basis_tools` command to convert the basis set. The conversion process preserved the phase/frequency convention fot he LCModel basis set. So to match that expected by FSL-MRS we also reverse it.\n",
" "
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%sx basis_tools convert example_data/example_mrsi/3T_slaser_32vespa_1250.BASIS 3T_slaser_32vespa_1250\n",
"%sx basis_tools conj 3T_slaser_32vespa_1250 3T_slaser_32vespa_1250"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Let's view the converted basis set"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"from fsl_mrs.utils import mrs_io\n",
"import matplotlib.pyplot as plt\n",
"basis = mrs_io.read_basis('3T_slaser_32vespa_1250')\n",
"fig = plt.figure(figsize=(8, 8))\n",
"_ = basis.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 2 - reference peak removal\n",
"LCModel basis sets have a reference peak in them at 0.0 ppm. This is removed by default in LCModel, but it isn't expected in FSL-MRS. Even the peak falls outside the default optimisation region, if you leave the reference peak in the summation across all metabolites can grow so large (as everything is scaled) that the peak distorts the baseline.\n",
"\n",
"We use the option `--remove_reference` to accomplish this. `--hlsvd` is added to make the removal more precise (at the cost of longer processing time)."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%sx basis_tools convert example_data/example_mrsi/3T_slaser_32vespa_1250.BASIS 3T_slaser_32vespa_1250_noTMS --remove_reference --hlsvd\n",
"%sx basis_tools conj 3T_slaser_32vespa_1250_noTMS 3T_slaser_32vespa_1250_noTMS\n",
"\n",
"basis = mrs_io.read_basis('3T_slaser_32vespa_1250_noTMS')\n",
"fig = plt.figure(figsize=(8, 8))\n",
"_ = basis.plot()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 3 - Remove duplicated basis spectra\n",
"No reference peaks any more! But we notice that there are two Gly basis sets. Are these actually identical?"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"import numpy as np\n",
"np.allclose(\n",
" basis.original_basis_array[:, basis.names.index('Gly_1')],\n",
" basis.original_basis_array[:, basis.names.index('Gly')])"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Ok, that's not good for the optimisation. We will remove the duplicate. Let's also say we aren't interested in 'HG' and we will also remove that metabolite."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%sx rm 3T_slaser_32vespa_1250_noTMS/Gly_1.json\n",
"%sx rm 3T_slaser_32vespa_1250_noTMS/HG.json\n",
"\n",
"basis = mrs_io.read_basis('3T_slaser_32vespa_1250_noTMS')\n",
"fig = plt.figure(figsize=(8, 8))\n",
"_ = basis.plot()\n",
"plt.xlim([5, 0])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Step 4 - Add default MM peaks\n",
"The above looks good now. But you notice that there aren't any macromolecular peaks, and this is a basis set for spectra with a 32 ms echo time. We therefore need to include suitable MM peaks. Ideally this is done with a measured (metabolite nulled) MM spectrum (see the example_usage/Example basis spectra creation.ipynb notebook), but we can use the FSL-MRS 'default' MM as a (poor) alternative.\n",
"\n",
"We do this by running the `basis_tools add_set` command with the `--add_MM` option.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": [
"%sx basis_tools add_set --add_MM 3T_slaser_32vespa_1250_noTMS 3T_slaser_32vespa_1250_noTMS_defaultMM\n",
"\n",
"basis = mrs_io.read_basis('3T_slaser_32vespa_1250_noTMS_defaultMM')\n",
"fig = plt.figure(figsize=(8, 8))\n",
"_ = basis.plot()\n",
"plt.xlim([5, 0])\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Conclusion\n",
"We've successfully converted an LCModel formatted basis set and added some default MM peaks. We can now use this in fitting of MRSI data!"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3.8.12 ('fsl_mrs')",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.8.12"
},
"orig_nbformat": 4,
"vscode": {
"interpreter": {
"hash": "6c11d37c82810953c5a08a185ec224dab920e965fab2a4fd7bf60d843c04e747"
}
}
},
"nbformat": 4,
"nbformat_minor": 2
}
%% Cell type:markdown id: tags:
# Command-line conversion of LCModel basis set
This notebook demonstrates the conversion of an existing LCModel basis set (`.basis` format).
In addition to the basic conversion step, three extra steps are performed:
1. Removal of the TMS peak from all spectra
2. Removal of duplicated and unwanted basis spectra
3. Addition of FSL-MRS default MM peaks
This converted spectrum will subsequently be used to fit the example MRSI data.
## Step 1 - conversion
We use the `basis_tools` command to convert the basis set. The conversion process preserved the phase/frequency convention fot he LCModel basis set. So to match that expected by FSL-MRS we also reverse it.
%% Cell type:code id: tags:
``` python
%sx basis_tools convert example_data/example_mrsi/3T_slaser_32vespa_1250.BASIS 3T_slaser_32vespa_1250
%sx basis_tools conj 3T_slaser_32vespa_1250 3T_slaser_32vespa_1250
```
%% Cell type:markdown id: tags:
Let's view the converted basis set
%% Cell type:code id: tags:
``` python
from fsl_mrs.utils import mrs_io
import matplotlib.pyplot as plt
basis = mrs_io.read_basis('3T_slaser_32vespa_1250')
fig = plt.figure(figsize=(8, 8))
_ = basis.plot()
```
%% Cell type:markdown id: tags:
## Step 2 - reference peak removal
LCModel basis sets have a reference peak in them at 0.0 ppm. This is removed by default in LCModel, but it isn't expected in FSL-MRS. Even the peak falls outside the default optimisation region, if you leave the reference peak in the summation across all metabolites can grow so large (as everything is scaled) that the peak distorts the baseline.
We use the option `--remove_reference` to accomplish this. `--hlsvd` is added to make the removal more precise (at the cost of longer processing time).
%% Cell type:code id: tags:
``` python
%sx basis_tools convert example_data/example_mrsi/3T_slaser_32vespa_1250.BASIS 3T_slaser_32vespa_1250_noTMS --remove_reference --hlsvd
%sx basis_tools conj 3T_slaser_32vespa_1250_noTMS 3T_slaser_32vespa_1250_noTMS
basis = mrs_io.read_basis('3T_slaser_32vespa_1250_noTMS')
fig = plt.figure(figsize=(8, 8))
_ = basis.plot()
```
%% Cell type:markdown id: tags:
## Step 3 - Remove duplicated basis spectra
No reference peaks any more! But we notice that there are two Gly basis sets. Are these actually identical?
%% Cell type:code id: tags:
``` python
import numpy as np
np.allclose(
basis.original_basis_array[:, basis.names.index('Gly_1')],
basis.original_basis_array[:, basis.names.index('Gly')])
```
%% Cell type:markdown id: tags:
Ok, that's not good for the optimisation. We will remove the duplicate. Let's also say we aren't interested in 'HG' and we will also remove that metabolite.
%% Cell type:code id: tags:
``` python
%sx rm 3T_slaser_32vespa_1250_noTMS/Gly_1.json
%sx rm 3T_slaser_32vespa_1250_noTMS/HG.json
basis = mrs_io.read_basis('3T_slaser_32vespa_1250_noTMS')
fig = plt.figure(figsize=(8, 8))
_ = basis.plot()
plt.xlim([5, 0])
plt.show()
```
%% Cell type:markdown id: tags:
## Step 4 - Add default MM peaks
The above looks good now. But you notice that there aren't any macromolecular peaks, and this is a basis set for spectra with a 32 ms echo time. We therefore need to include suitable MM peaks. Ideally this is done with a measured (metabolite nulled) MM spectrum (see the example_usage/Example basis spectra creation.ipynb notebook), but we can use the FSL-MRS 'default' MM as a (poor) alternative.
We do this by running the `basis_tools add_set` command with the `--add_MM` option.
%% Cell type:code id: tags:
``` python
%sx basis_tools add_set --add_MM 3T_slaser_32vespa_1250_noTMS 3T_slaser_32vespa_1250_noTMS_defaultMM
basis = mrs_io.read_basis('3T_slaser_32vespa_1250_noTMS_defaultMM')
fig = plt.figure(figsize=(8, 8))
_ = basis.plot()
plt.xlim([5, 0])
plt.show()
```
%% Cell type:markdown id: tags:
## Conclusion
We've successfully converted an LCModel formatted basis set and added some default MM peaks. We can now use this in fitting of MRSI data!
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
This source diff could not be displayed because it is stored in LFS. You can view the blob instead.
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment