Analysis Primer

From VISTA LAB WIKI

(Difference between revisions)
Jump to: navigation, search
Line 1: Line 1:
-
There are some methods (e.g., coordinate frame alignment) that are common to all of the analysis techniques. We describe these techniques here.  
+
There are some methods (e.g., coordinate frame alignment) that are common to all of the analysis techniques. We describe these techniques here. There are some methods (e.g., coordinate frame alignment) that are common to all of the analysis techniques. We describe these techniques here.  
-
There are some methods (e.g., coordinate frame alignment) that are common to all of the analysis techniques. We describe these techniques here.
+
-
== Create Time Series from k-space data ==
+
== Create Time Series from k-space data ==
-
Functional measurements consist of time series of responses. The first step in analyzing the functional data is to create this time series.
+
Functional measurements consist of time series of responses. The first step in analyzing the functional data is to create this time series.  
For stock vendor sequences (such as EPI), it is fully automatic and the resulting images appear in the vendor's image database. You can get these images from that database as DICOM files at the scanner.  
For stock vendor sequences (such as EPI), it is fully automatic and the resulting images appear in the vendor's image database. You can get these images from that database as DICOM files at the scanner.  
-
For non-vendor sequences, such as the spiral methods favored at the Lucas Center, you may need to run a separate 'recon' program to reconstruct the images from the raw k-space data that the scanner collects. At Stanford's Lucas Center, we often use Gary Glover's spiral sequence, which by default runs Gary's 'grecons' program to convert the raw k-space data (GE's 'P-file') and produce:
+
For non-vendor sequences, such as the spiral methods favored at the Lucas Center, you may need to run a separate 'recon' program to reconstruct the images from the raw k-space data that the scanner collects. At Stanford's Lucas Center, we often use Gary Glover's spiral sequence, which by default runs Gary's 'grecons' program to convert the raw k-space data (GE's 'P-file') and produce:  
-
# PNNNNN.7.mag: raw magnitude image data
+
#PNNNNN.7.mag: raw magnitude image data  
-
# PNNNNN.7.hdr: dump of raw GE header info
+
#PNNNNN.7.hdr: dump of raw GE header info  
-
# E...PNNNNN.7: human-readable header info
+
#E...PNNNNN.7: human-readable header info
-
These data files do not go into the scanner's image database; instead they are transferred to the analysis computer by file transfer protocols. At the Lucas center, we transfer these files directly from the scanner's 'raw' data directory.
+
These data files do not go into the scanner's image database; instead they are transferred to the analysis computer by file transfer protocols. At the Lucas center, we transfer these files directly from the scanner's 'raw' data directory.  
-
== Reformat raw images ==
+
== Reformat raw images ==
-
See: '''mrLoad''' and '''mrInit''' in [[Initialization]]
+
See: '''mrLoad''' and '''mrInit''' in [[Initialization]]  
-
The functional MRI program [[Initialization#mrInit|mrInit]] reads E and P.mag files and generates a set of Matlab, (.mat) files. We intend to move to the more standard [http://nifti.nimh.nih.gov/nifti-1 NIFTI] file format. The fMRI timeseries for each scan will then be stored in one 4-d NIFTI file. The code to manage different file types has been written (Rory) and is fairly advanced. We plan to integrate this function, mrLoad, from mrVista2.
+
The functional MRI program [[Initialization#mrInit|mrInit]] reads E and P.mag files and generates a set of Matlab, (.mat) files. We intend to move to the more standard [http://nifti.nimh.nih.gov/nifti-1 NIFTI] file format. The fMRI timeseries for each scan will then be stored in one 4-d NIFTI file. The code to manage different file types has been written (Rory) and is fairly advanced. We plan to integrate this function, mrLoad, from mrVista2.  
-
The dMRI analysis pipeline (mrDiffusion) expects a 4-d NIFTI file of the raw diffusion-weighted (DW) images, including a correct scanner-space transform encoded in the header. To process dMRI data, you also need to know the DW strength (specified as the 'b-value') and the list of DW directions that correspond to the raw DW images. mrDiffusion expects these to be stored in b-vals and b-vecs files, the same format that [http://www.fmrib.ox.ac.uk/fsl/fdt/fdt_intro.html FSL's fdt] uses.
+
The dMRI analysis pipeline (mrDiffusion) expects a 4-d NIFTI file of the raw diffusion-weighted (DW) images, including a correct scanner-space transform encoded in the header. To process dMRI data, you also need to know the DW strength (specified as the 'b-value') and the list of DW directions that correspond to the raw DW images. mrDiffusion expects these to be stored in b-vals and b-vecs files, the same format that [http://www.fmrib.ox.ac.uk/fsl/fdt/fdt_intro.html FSL's fdt] uses.  
-
== Slice timing correction ==
+
== Slice timing correction ==
-
See: '''AdjustSliceTiming([scans], [typeName])''' and  
+
See: '''AdjustSliceTiming([scans], [typeName])''' and '''mrSliceTiming(ts,scan,slice,method)'''  
-
'''mrSliceTiming(ts,scan,slice,method)'''
+
Slice time correction corrects for differences in the slice-specific acquisition time.  
Slice time correction corrects for differences in the slice-specific acquisition time.  
-
Slices are acquired generally either in one of three orderings:
+
Slices are acquired generally either in one of three orderings:  
-
* ascending: [1 2 3 4 5 6 7 8 9 10]
+
*ascending: [1 2 3 4 5 6 7 8 9 10]  
-
* descending: [10 9 8 7 6 5 4 3 2 1]
+
*descending: [10 9 8 7 6 5 4 3 2 1]  
-
* interleaved: [1 3 5 7 9 2 4 6 8 10]
+
*interleaved: [1 3 5 7 9 2 4 6 8 10]
-
If you acquire a volume of data every 1.5 seconds (TR=1.5s) the delay between the first and last slice is the TR. The delay between adjacent slices is either TR/nSlices or TR/2, depending on whether you use an ordered or interleaved method.  
+
If you acquire a volume of data every 1.5 seconds (TR=1.5s) the delay between the first and last slice is the TR. The delay between adjacent slices is either TR/nSlices or TR/2, depending on whether you use an ordered or interleaved method.  
These differences are modest compared to the hemodynamic response function (4-6 sec). For many purposes, it is not essential to correct the slice timing.  
These differences are modest compared to the hemodynamic response function (4-6 sec). For many purposes, it is not essential to correct the slice timing.  
Line 43: Line 41:
It is possible to approach slice timing correction two ways.  
It is possible to approach slice timing correction two ways.  
-
First, you can interpolate the raw data to correct for these timing differences. The advantage of this approach is that later preprocessing steps, such as smoothing or interpolating across slices (e.g. due to motion correction, spatial normalization) are not affected by slice-time differences. The disadvantage is that the data need to be interpolated and extrapolated at the edges. The slice timing also smooths the data a slight amount by the interpolation.  
+
First, you can interpolate the raw data to correct for these timing differences. The advantage of this approach is that later preprocessing steps, such as smoothing or interpolating across slices (e.g. due to motion correction, spatial normalization) are not affected by slice-time differences. The disadvantage is that the data need to be interpolated and extrapolated at the edges. The slice timing also smooths the data a slight amount by the interpolation.  
-
You will need to chose an interpolation method. At the Lucas Center (a) the slice order is ascending and (b) the slice acquisition times are equally spaced in time. These assumptions are not necessarily valid for other sequences. You can check the slice timing sequence for a data set by XXX???XXX. This will be incorporated properly at some point - ask BW/Bob
+
You will need to chose an interpolation method. At the Lucas Center (a) the slice order is ascending and (b) the slice acquisition times are equally spaced in time. These assumptions are not necessarily valid for other sequences. You can check the slice timing sequence for a data set by XXX???XXX. This will be incorporated properly at some point - ask BW/Bob  
-
Second, you can correct for slice timing differences in the final analysis. In a GLM analysis your design matrix can be adjusted for each slice and in a phase-encoded design the phase of the best fitting sinusoid can be adjusted according to the acquisition of each slice. This has the advantage that no interpolation of the raw data is required. The disadvantage is that this correction is late in the processing pipeline and the slice timing may be altered by smoothing and spatial resampling of the data (e.g. due to motion correction or spatial normalization).
+
Second, you can correct for slice timing differences in the final analysis. In a GLM analysis your design matrix can be adjusted for each slice and in a phase-encoded design the phase of the best fitting sinusoid can be adjusted according to the acquisition of each slice. This has the advantage that no interpolation of the raw data is required. The disadvantage is that this correction is late in the processing pipeline and the slice timing may be altered by smoothing and spatial resampling of the data (e.g. due to motion correction or spatial normalization). (I don't think we ever do this. This comment should be removed, I think -- BW).  
-
(I don't think we ever do this. This comment should be removed, I think -- BW).
+
-
When we use multiple shots, the ordering is [1a 2a 3a ... 1b 2b 3b ...] where a and b are two shots. So the effective spacing between the acquisition times for slice timing = [1 / number of slices / number of shots]. If you have participants who are likely to move (e.g., children), you are probably better off using a single shot. However, you may be working with data from other labs that use multiple shots, or you may have subjects that are not likely to move for which multiple shots would result in nicer-looking data (e.g., Logothetis' 8-shot scans of anaesthetized monkeys).
+
When we use multiple shots, the ordering is [1a 2a 3a ... 1b 2b 3b ...] where a and b are two shots. So the effective spacing between the acquisition times for slice timing = [1 / number of slices / number of shots]. If you have participants who are likely to move (e.g., children), you are probably better off using a single shot. However, you may be working with data from other labs that use multiple shots, or you may have subjects that are not likely to move for which multiple shots would result in nicer-looking data (e.g., Logothetis' 8-shot scans of anaesthetized monkeys).  
-
== Smoothing ==
+
== Smoothing ==
-
To smooth or not to smooth...
+
To smooth or not to smooth...  
-
=== Spatial smoothing ===
+
=== Spatial smoothing ===
-
=== Temporal smoothing ===
+
=== Temporal smoothing ===
-
== Time series motion correction ==
+
== Time series motion correction ==
Many labs do this routinely on every dataset. In the VISTA lab, we usually don't do this unless needed. There are several methods implemented in mrVista.  
Many labs do this routinely on every dataset. In the VISTA lab, we usually don't do this unless needed. There are several methods implemented in mrVista.  
-
For dMRI, you will probably need to do eddy-current correction, which might also include motion correction. In fact, the same exact algorithms can be use for fMRI motion correction and dMRI eddy/motion correction. The only differences are that you migth want to limit fMRI motion correction to a rigid-body (6-parameter) transform, while eddy current distortions require more degrees of freedom (e.g., an affine (12-parameter) transform or specialized constrained non-linear transforms). Also, you must use a robust error estimator, such as mutual information, since the dMRI image contrast varies with the different diffusion weighing strengths and directions. In mrDiffusion, we use our implementation of the Rohde et. al. (2004 in MRM) 14-parameter transform to do motion and eddy-current correction (see the section on [[DTI#dtiRaw_pre-processing_pipeline|the mrDiffusion pre-processing pipeline]]).
+
For dMRI, you will probably need to do eddy-current correction, which might also include motion correction. In fact, the same exact algorithms can be use for fMRI motion correction and dMRI eddy/motion correction. The only differences are that you migth want to limit fMRI motion correction to a rigid-body (6-parameter) transform, while eddy current distortions require more degrees of freedom (e.g., an affine (12-parameter) transform or specialized constrained non-linear transforms). Also, you must use a robust error estimator, such as mutual information, since the dMRI image contrast varies with the different diffusion weighing strengths and directions. In mrDiffusion, we use our implementation of the Rohde et. al. (2004 in MRM) 14-parameter transform to do motion and eddy-current correction (see the section on [[DTI#dtiRaw_pre-processing_pipeline|the mrDiffusion pre-processing pipeline]]).  
-
The basic steps for motion correction are:
+
The basic steps for motion correction are:  
-
*Estimate motion
+
 
-
*Resample images to correct the motion
+
*Estimate motion  
 +
*Resample images to correct the motion  
*Find and remove bad images?
*Find and remove bad images?
-
== fMRI slow-drift trend removal ==
+
== fMRI slow-drift trend removal ==
-
For most fMRI sequences, the MR image intensity will usually drift over time. This is often attributed to thermal drift as the gradient heat up, but may have other causes as well. (Drift in dMRI sequences?)
+
For most fMRI sequences, the MR image intensity will usually drift over time. This is often attributed to thermal drift as the gradient heat up, but may have other causes as well. (Drift in dMRI sequences?)  
-
One approach to dealing with such global signal drift is to simply add a slow-drift term to the model that you fit to the data (see below). In mrVista
+
One approach to dealing with such global signal drift is to simply add a slow-drift term to the model that you fit to the data (see below). In mrVista we often estimate the drift separate from the model fit and explicitly remove it from the time series.  
-
we often estimate the drift separate from the model fit and explicitly remove it from the time series.
+
-
== Compute within-subject alignment ==
+
== Compute within-subject alignment ==
[[Image:HeaderAlign.png|thumb|left|300px|An EPI image (middle) roughly aligned with a T1-weighted high-res anatomy of the same brain (top) based on scanner header information.]]
[[Image:HeaderAlign.png|thumb|left|300px|An EPI image (middle) roughly aligned with a T1-weighted high-res anatomy of the same brain (top) based on scanner header information.]]
[[Image:AutoAlign.png|thumb|left|300px|Alignment refined by a mutual information algorithm. In the bottom row, the EPI is shown as a red overlay on T1-weighted image.]]
[[Image:AutoAlign.png|thumb|left|300px|Alignment refined by a mutual information algorithm. In the bottom row, the EPI is shown as a red overlay on T1-weighted image.]]
-
The high-resolution anatomical scan (usually a T1-weighted SPGR or MP-RAGE) forms a common reference space for all the data from an individual subject. By using the reference space, we can combine data across different fMRI scan sessions and between fMRI and dMRI datasets.
+
The high-resolution anatomical scan (usually a T1-weighted SPGR or MP-RAGE) forms a common reference space for all the data from an individual subject. By using the reference space, we can combine data across different fMRI scan sessions and between fMRI and dMRI datasets.  
 +
 
 +
<br>
 +
 
-
<p>
+
== Compute alignment to a standard space ==
-
== Compute alignment to a standard space ==
+
In the VISTA lab, we routinely ac-pc align each subject's high-res anatomical. This provides a rough alignment across subjects, useful for finding similar brain regions. However, we usually combine data across subjects using ROI analysis methods.  
In the VISTA lab, we routinely ac-pc align each subject's high-res anatomical. This provides a rough alignment across subjects, useful for finding similar brain regions. However, we usually combine data across subjects using ROI analysis methods.  
-
* ROI analysis
+
*ROI analysis  
-
* alignment to a template (e.g., MNI space)
+
*alignment to a template (e.g., MNI space)
-
== Fit a model to the data ==
+
== Fit a model to the data ==
To interpret the data, you will need to fit a model.  
To interpret the data, you will need to fit a model.  
-
For fMRI data, many possibilities exist. For example, you might fit a simple sinusoidal model for classic retinotopy data and simple block designs. This model produces three parameters per voxel- the sinusoid amplitude and phase, and the coherence- which is the amplitude normalized by the other frequency componenets and is comparable to a correlation coefficient.  
+
For fMRI data, many possibilities exist. For example, you might fit a simple sinusoidal model for classic retinotopy data and simple block designs. This model produces three parameters per voxel- the sinusoid amplitude and phase, and the coherence- which is the amplitude normalized by the other frequency componenets and is comparable to a correlation coefficient.  
For more complicated block-design paradigms and event-related designs, you may build a general linear model (GLM) that is very similar to that used in SPM.  
For more complicated block-design paradigms and event-related designs, you may build a general linear model (GLM) that is very similar to that used in SPM.  
-
Finally, Serge Dumoulin has developed a more sophisticated model for retinotopy analysis- the 'pRF' model.
+
Finally, Serge Dumoulin has developed a more sophisticated model for retinotopy analysis- the 'pRF' model.  
-
For diffusion data (dMRI), the most common model is the diffusion-tensor (DTI). mrDiffusion currently fits the diffusion tensor using a simple least-sqaures approach; we are exploring more robust tensor-fitting methods.
+
For diffusion data (dMRI), the most common model is the diffusion-tensor (DTI). mrDiffusion currently fits the diffusion tensor using a simple least-sqaures approach; we are exploring more robust tensor-fitting methods.  
-
== Infer something interesting ==
+
== Infer something interesting ==
Finally, you'll look at your data and the resulting model parameters, combine data across some subjects, test hypotheses, make a discovery, and publish ground-breaking papers.
Finally, you'll look at your data and the resulting model parameters, combine data across some subjects, test hypotheses, make a discovery, and publish ground-breaking papers.

Revision as of 16:14, 17 July 2015

There are some methods (e.g., coordinate frame alignment) that are common to all of the analysis techniques. We describe these techniques here. There are some methods (e.g., coordinate frame alignment) that are common to all of the analysis techniques. We describe these techniques here.

Contents

Create Time Series from k-space data

Functional measurements consist of time series of responses. The first step in analyzing the functional data is to create this time series.

For stock vendor sequences (such as EPI), it is fully automatic and the resulting images appear in the vendor's image database. You can get these images from that database as DICOM files at the scanner.

For non-vendor sequences, such as the spiral methods favored at the Lucas Center, you may need to run a separate 'recon' program to reconstruct the images from the raw k-space data that the scanner collects. At Stanford's Lucas Center, we often use Gary Glover's spiral sequence, which by default runs Gary's 'grecons' program to convert the raw k-space data (GE's 'P-file') and produce:

  1. PNNNNN.7.mag: raw magnitude image data
  2. PNNNNN.7.hdr: dump of raw GE header info
  3. E...PNNNNN.7: human-readable header info

These data files do not go into the scanner's image database; instead they are transferred to the analysis computer by file transfer protocols. At the Lucas center, we transfer these files directly from the scanner's 'raw' data directory.

Reformat raw images

See: mrLoad and mrInit in Initialization

The functional MRI program mrInit reads E and P.mag files and generates a set of Matlab, (.mat) files. We intend to move to the more standard NIFTI file format. The fMRI timeseries for each scan will then be stored in one 4-d NIFTI file. The code to manage different file types has been written (Rory) and is fairly advanced. We plan to integrate this function, mrLoad, from mrVista2.

The dMRI analysis pipeline (mrDiffusion) expects a 4-d NIFTI file of the raw diffusion-weighted (DW) images, including a correct scanner-space transform encoded in the header. To process dMRI data, you also need to know the DW strength (specified as the 'b-value') and the list of DW directions that correspond to the raw DW images. mrDiffusion expects these to be stored in b-vals and b-vecs files, the same format that FSL's fdt uses.

Slice timing correction

See: AdjustSliceTiming([scans], [typeName]) and mrSliceTiming(ts,scan,slice,method)

Slice time correction corrects for differences in the slice-specific acquisition time.

Slices are acquired generally either in one of three orderings:

  • ascending: [1 2 3 4 5 6 7 8 9 10]
  • descending: [10 9 8 7 6 5 4 3 2 1]
  • interleaved: [1 3 5 7 9 2 4 6 8 10]

If you acquire a volume of data every 1.5 seconds (TR=1.5s) the delay between the first and last slice is the TR. The delay between adjacent slices is either TR/nSlices or TR/2, depending on whether you use an ordered or interleaved method.

These differences are modest compared to the hemodynamic response function (4-6 sec). For many purposes, it is not essential to correct the slice timing.

It is possible to approach slice timing correction two ways.

First, you can interpolate the raw data to correct for these timing differences. The advantage of this approach is that later preprocessing steps, such as smoothing or interpolating across slices (e.g. due to motion correction, spatial normalization) are not affected by slice-time differences. The disadvantage is that the data need to be interpolated and extrapolated at the edges. The slice timing also smooths the data a slight amount by the interpolation.

You will need to chose an interpolation method. At the Lucas Center (a) the slice order is ascending and (b) the slice acquisition times are equally spaced in time. These assumptions are not necessarily valid for other sequences. You can check the slice timing sequence for a data set by XXX???XXX. This will be incorporated properly at some point - ask BW/Bob

Second, you can correct for slice timing differences in the final analysis. In a GLM analysis your design matrix can be adjusted for each slice and in a phase-encoded design the phase of the best fitting sinusoid can be adjusted according to the acquisition of each slice. This has the advantage that no interpolation of the raw data is required. The disadvantage is that this correction is late in the processing pipeline and the slice timing may be altered by smoothing and spatial resampling of the data (e.g. due to motion correction or spatial normalization). (I don't think we ever do this. This comment should be removed, I think -- BW).

When we use multiple shots, the ordering is [1a 2a 3a ... 1b 2b 3b ...] where a and b are two shots. So the effective spacing between the acquisition times for slice timing = [1 / number of slices / number of shots]. If you have participants who are likely to move (e.g., children), you are probably better off using a single shot. However, you may be working with data from other labs that use multiple shots, or you may have subjects that are not likely to move for which multiple shots would result in nicer-looking data (e.g., Logothetis' 8-shot scans of anaesthetized monkeys).

Smoothing

To smooth or not to smooth...

Spatial smoothing

Temporal smoothing

Time series motion correction

Many labs do this routinely on every dataset. In the VISTA lab, we usually don't do this unless needed. There are several methods implemented in mrVista.

For dMRI, you will probably need to do eddy-current correction, which might also include motion correction. In fact, the same exact algorithms can be use for fMRI motion correction and dMRI eddy/motion correction. The only differences are that you migth want to limit fMRI motion correction to a rigid-body (6-parameter) transform, while eddy current distortions require more degrees of freedom (e.g., an affine (12-parameter) transform or specialized constrained non-linear transforms). Also, you must use a robust error estimator, such as mutual information, since the dMRI image contrast varies with the different diffusion weighing strengths and directions. In mrDiffusion, we use our implementation of the Rohde et. al. (2004 in MRM) 14-parameter transform to do motion and eddy-current correction (see the section on the mrDiffusion pre-processing pipeline).

The basic steps for motion correction are:

  • Estimate motion
  • Resample images to correct the motion
  • Find and remove bad images?

fMRI slow-drift trend removal

For most fMRI sequences, the MR image intensity will usually drift over time. This is often attributed to thermal drift as the gradient heat up, but may have other causes as well. (Drift in dMRI sequences?)

One approach to dealing with such global signal drift is to simply add a slow-drift term to the model that you fit to the data (see below). In mrVista we often estimate the drift separate from the model fit and explicitly remove it from the time series.

Compute within-subject alignment

An EPI image (middle) roughly aligned with a T1-weighted high-res anatomy of the same brain (top) based on scanner header information.
Alignment refined by a mutual information algorithm. In the bottom row, the EPI is shown as a red overlay on T1-weighted image.

The high-resolution anatomical scan (usually a T1-weighted SPGR or MP-RAGE) forms a common reference space for all the data from an individual subject. By using the reference space, we can combine data across different fMRI scan sessions and between fMRI and dMRI datasets.



Compute alignment to a standard space

In the VISTA lab, we routinely ac-pc align each subject's high-res anatomical. This provides a rough alignment across subjects, useful for finding similar brain regions. However, we usually combine data across subjects using ROI analysis methods.

  • ROI analysis
  • alignment to a template (e.g., MNI space)

Fit a model to the data

To interpret the data, you will need to fit a model.

For fMRI data, many possibilities exist. For example, you might fit a simple sinusoidal model for classic retinotopy data and simple block designs. This model produces three parameters per voxel- the sinusoid amplitude and phase, and the coherence- which is the amplitude normalized by the other frequency componenets and is comparable to a correlation coefficient.

For more complicated block-design paradigms and event-related designs, you may build a general linear model (GLM) that is very similar to that used in SPM.

Finally, Serge Dumoulin has developed a more sophisticated model for retinotopy analysis- the 'pRF' model.

For diffusion data (dMRI), the most common model is the diffusion-tensor (DTI). mrDiffusion currently fits the diffusion tensor using a simple least-sqaures approach; we are exploring more robust tensor-fitting methods.

Infer something interesting

Finally, you'll look at your data and the resulting model parameters, combine data across some subjects, test hypotheses, make a discovery, and publish ground-breaking papers.

Personal tools