Home
mrVista WorkFLow
Scanning
Post Scanning
fMRI Data Analysis
Original Lab Manual
Other Resources









 

fMRI Data Analysis

mrLoadRet-3.0 User Guide

Overview:

mrLoadRet is our in-house fMRI data analysis program. It offers 3 views of the cortex: Inplane (from each functional scanning session), Volume or Gray (the gray view restricts the data to the segmented gray matter), and Flat (flattened gray matter). This combination of views allows several functions: viewing functional data as 'maps' on any of these views of cortex; defining regions-of-interest (ROIs); plotting responses for ROIs; and other statistical/data analysis.
mrLoadRet has a "plug-in": mrAlign, which allows either automatic or manual alignment of the inplanes to the reference volume.

mrLoadRet is not any of the following: mrGray (for segmenting gray matter), mrFlatMesh (for unfolding gray matter), mrInitRet (for initial data reconstruction). These are separate applications; see the corresponding labmanual pages for more information.


Directory structure:

The top-level data directory is typically named by date (e.g., 101998). This directory contains a couple critical files:
Readme or Readme.txt:
mrSESSION.mat: contains the mrSESSION struction that is created by mrInitRet and that is loaded by mrLoadRet. This structure has fields for all the information that mrLoadRet needs to handle this data set. See mrLoadRet hacking for details.
AlignParams.mat, bestrotvol.mat: These files contain information about the alignment between the data from this session and the subject's volume anatomy.
There are subdirectories for each view (Inplane, Flat, Gray, and Volume). There may be multiple Flat subdirectories, each corresponding to a different unfold, but they must all be called 'Flat-*'. There is also a Raw subdirectory that contains the raw Pfiles, inplane anatomy images and screen save.
Within each of the view subdirectories, there is an ROIs subdirectory containing various ROI files (e.g., V1.mat).

Within each of the view subdirectories, there are one or more data subdirectories for each of the data types. If there is only one data subdirectory, it should be called 'Original'. Other data types might include 'Averages' (that is created by averaging data across 2 or more scans within the session) or 'Complex' (for complex-valued recons).

Within each of the data subdirectories, there will typically be a corAnal.mat file along with other maps (e.g., meanMap.mat). and a TSeries directory for the time-series files generated by the recon process. In TSeries, there is a separate subdir for each scan (Scan1, Scan2, etc), and in each of these, there is a tSeries*.mat file, 1 for each slice. So the full path to a particular tSeries file will be something like: <date>/Inplane/Original/TSeries/Scan1/tSeries1.mat

There may also be an Analysis subdirectory within one of the data subdirectories which, if it exists, contains the results of some data analysis procedure applied to this session's data (see Analysis of Blocked Alternation Data Sets below).

Getting Started:

Type "mrLoadRet" at the MATLAB prompt. The mrLoadRet Inplane window will open. If mrLoadRet doesn't start up right, check what directory you are in. It should be the top-level directory that includes mrSESSION.
From the File menu, select "Load Anatomies" and "Load Correlation Matrices". Note that the sliders on the top select the scan and slice.
From the View menu, you can now look at a map of your functional data overlaid on the inplanes. You can view different aspects fo the data: by correlation, by amplitude, by phase, etc.
Cothresh and phWindow. The 3 sliders on the right of the window are labeled cothresh and phWindow min/max. Adjust them to reasonable levels: standard practice is a cothresh of ~.4, and a phWindow of 0-pi/2 (assuming that you're looking for in-phase responses). Setting these sliders will remove the functional data map points that do not satisfy these correlation/phase criteria. Sliding any slider to its maximum value allows you to input the desired value at the MATLAB prompt.
ROIs can be created and modified using tools in the ROI menu. Pre-existing ROIs can also be loaded using the File menu. You can then restrict the ROI according to the current settings of the cothresh and phWindow sliders. If you have more than one ROI defined, use the ROI popup to choose which ROI is currently selected. Note that several keyboard shortcuts ('Ctrl-key') have been added to the ROIs menu, and to other frequently-used commands. They are noted in the menu items.
Plots: From the Plot/Across Scans menu, select one of the options. Note that for the projected amplitude plots, the current scan is used to compute the reference ohase vector. From the Plot/Selected Scan menu, select one of the options. Resulting plot shows data averaged across pixels within the current ROI.
Use the data type popup to choose the data type (e.g., 'Original' or 'Averages'). After switching data types, you have to reload the data because the correlation (and other parameter) maps are different for each data type.
From the Analysis menu, you can compute the correlation, amplitude, and phase maps along with several other pre-defined parameter maps.
From the Window menu, you can open another Inplane view. Each view may be set to the same or a different data type.
From the Color Map menu, you can select different color maps for the pseudo-color representations of the data.

Advanced Features:

Selecting a view. Because mrLoadRet-3 now allows you have multiple Inplane windows open simultaneoulsy, one of them is designated as the selected Inplane (likewise for Volume/Gray and Flat windows). This is critical, for example, when transforming data from the Inplane to the Gray view. The data are grabbed from the selected Inplane, not from any of the other inplanes. To specify a window as the selected Inplane, just click on it almost anywhere within the window (clicking on the title bar won't do it), and likewise for the Volume/Gray and Flat windows.
Opening Volume, Gray, and Flat views. From the Window menu, you can open one or more Volume or Gray views, and/or one or more Flat views. Each Flat view may correspond to a different unfold. But all Gray and Flat views must correspond to the same gray matter segmentation and the same volume anatomy (vAnatomy.mat). The first time you open a Volume or Gray view, mrLoadRet will look for the vAnatomy file in a standard location (for windows, x:\anatomy\<subject>; for unix/linux, /usr/local/mri/anatomy/<subject>). If no vAnatomy is found, mrLoadRet will prompt you to locate the file. The first time you open a Gray view, mrLoadRet will prompt you for the gray classification file (i.e., the output generated by mrGray). When you open a Flat view (from the Window menu), if there is more than Flat-* subdirectory, mrLoadRet will prompt you to pick one. The first time you open a Flat view, mrLoadRet will prompt you to locate the flat.mat file (i.e., the output of mrUnfold).

From the segmentation menu you can:

Re-install the segmentation. This deletes all of the data files in all of the Gray and Flat subdirectories (except for the ROIs), and prompts you to locate the gray classification file. This is useful when you redo the segmentation in mrGray or if, for some reason, the Gray and Flat views get out of sync.
Re-install an existing unfold. This deletes all of the data files in the current Flat subdirectory (except optionally the ROIs), and prompts you to locate the flat.mat file.
Install a new unfold. Makes a new Flat-* subdirectory and prompts you for the new flat.mat file.
Get information about the segmentation and unfold and check to make sure that the current flat view is consistent with (i.e., was unfolded from) the right segmentation.
Transforming data across views. From the Xform menu you can transform data (parameter maps like co, amp, ph, meanMap) from the Inplane view to the Volume or Gray views and then onto the Flat view. You can also transform tSeries from Inplane to Gray (but not to Volume because the data files would get too big) and then onto the Flat. Because the data start off in the Inplane, it does not make sense to transform data in the other direction (e.g., from Flat to Gray to Inplane). You can transform ROIs in either both directions (Inplane <-> Volume/Gray <-> Flat) because some ROIs are defined on the Flat map. The xform menu on any View lets you do operations for which that view is the destination. For example, to transform an ROI from the Inplane to the Volume, you have to have an Inplane window open and selected (see Selecting a view above) with a currently selected ROI. The Xform/Current ROI will tranform the currently selected ROI from the Inplane view to the Volume view.
The Gray view should be considered home base. This is the place where data and ROIs from multiple scanning sessions can easily be combined or analyzed concurrently. I recommend that you do most/all of your data analysis in this view. The data in the Gray view are identical to those in the Inplance view (just oversampled because of the smaller voxel size in the Volume/Gray). The ROIs are slightly different in the two views because they get "fat" when transformed from the Gray to the Inplane (because the Inplane voxels are bigger).

The Flat view should be considered a visualization tool. Most forms of data analysis should NOT be performed in the Flat view because the mapping from Gray to Flat is neither 1-to-1 nor onto. There are many Flat pixels that collapse across several gray nodes (e.g., the ~4mm thickness of the gray matter will collapse to a single flat pixel). There are many Flat pixels that don't contain any gray nodes, in part because of distortion inherent in the unfolding.

Data Types. mrLoadRet-3.0 supports multiple data types. Typically, a session will include the 'Original' data type. Other data types might include 'Averages' (that is created by averaging data across 2 or more scans within the session) or 'Complex' (for complex-valued recons) or 'Copied' (for scans that are copied in from other scanning sessions). To switch between data types, use the pull down menu in the top-right corner of each window. This clears all of the data fields so, for example, you need to load the corAnal matrices after switching data types. To transform data from one view to another, the two views must be set to the same data type.

Motion correction. There are several ways that you can check for motion artifacts or other artifacts in the images. First, view the tSeries movie (from the View menu) to see if the subject moved during a functional scan. A second way to check for motion during a scan is to plot either/both the rmse across frames and/or the max difference across frames (from the Plots menu). Third, compute the mean map (from the Analysis menu) and view it overlaid (with appropriate slider settings) superimposed on the inplane anatomies to make sure that the subject didn't move in between scans.You can apply motion correction if necessary (from the Analysis menu).

Checking for other artifacts. The tSeries Pict option from the Analysis menu makes a picture of the (temporal) average of the functional images. Certain image artifacts are easily visible in these pictures (which should look basically like brain slices). Gary likes to look at maps of the stdev over time (Analysis/compute stdMap) to check for artifacts, which should look like a random field. Some artifacts show up as swirly patterns or smooth gradients in the phase maps, across the whole slice. It is generally a good idea to check for image artifacts in each and every data set.

Copying/Merging data from multiple scanning sessions. The function copyScan can be used copy the tSeries data from one scanning session to another. This is possible only by working through the Gray view (shared coordinates for both scanning sessions). Here's how to do it:

Run mrAlign on both scanning sessions
Run mrLoadRet in the session1, the session that you want to copy from
Open a Gray window (Window/Open gray window). Note that segmentation must have already been done.
Copy one or more tSeries from Inplane to Gray (Xform/Inplane->Volume/tseries)
Quit mrLoadRet (close all; clear all)
Run mrLoadRet in the session2, the session that you want to copy to
Open a Gray window.
Call the copyScan function by typing at the matlab prompt (sorry but there is no menu interface to this as of now):
copyScan(session1Path, dataType, scanNumber);
where
session1Path is a string specifying the full path to session1, e.g., 'e:\mri\test'
dataType is the data type of the scan you want to copy, e.g., 'Original' or 'Averages'
scanNumber is the scan number, e.g., 1 or 2

The copied scan shows up under the "Copied" dataType.

You can also make a new "virtual" scanning session and copy in data from several other scanning sessions, using the functions mergeSessions and copyScan. Here's an example script illustrating how to do this:

sessionDirs = {'e:/mri/test', 'j:/djhRetinotopy/030701a'};
newSession = 'e:/mri/foo';
mergeSessions(sessionDirs,newSession);
chdir(newSession)

mrLoadRet
copyScan(sessionDirs{1},'Original',1);
copyScan(sessionDirs{1},'Original',2);
copyScan(sessionDirs{2},'Averages',1);
copyScan(sessionDirs{2},'Averages',2);
The new session directory (e:/mri/foo in the example) must already exist before executing this script. See below for more examples of how to write scripts.

Analysis:

Standard correlation analysis. The function computeCorAnal (called when you choose compute corAnal from the Analysis menu) does our standard data analysis for periodic stimulation protocols, including retintopy measurements. It "fits" a sinusoid to the time-series at each pixel and returns 3 numbers: the amplitude of the best-fit sinusoid, the phase of the best-fit sinusoid, and the correlation between the time-series and the best-fit sinusoid. The amplitude is a measure of the signal strength, the phase reflects the hemodynamic delay, and the correlation is a measure of signal-to-noise (see any number of our papers for details).
You can compute the corAnal in the Inplane view and transform the results to the Gray view, or you can transform the tSeries to the Gray view first and then compute the corAnal. Either way will produce the same results.

Detrending and inhomogeneity correction. In analyzing our data, we typically remove the slow drift in the fMRI signal. We also typically divide by the mean to convert to units of percent BOLD signal change (rather than the raw intensity units, and to compensate for the distance from the surface coil. [The function loadtSeries returns the raw tSeries. The function percentTSeries does the detrending and inhomogeneity correction and holds the result in the view.tSeries field.]

There are several options for how to remove the baseline, depending on the value of mrSESSION.analysis.detrend:

0 no trend removal
1 highpass trend removal (default)
2 quartic removal
-1 linear trend removal
There are two options for how to compensate for distance from the coil, depending on the value of mrSESSION.analysis.inhomoCorrection:
0 do nothing
1 divide by the mean, independently at each voxel (default)
2 divide by the null condition (not yet implemented)
3 divide by anything you like, e.g., a robust estimate of intensity inhomogeneity
The different inhomogeneity correction options give identical correlation and phase values. The difference is evident only in the amplitudes. For inhomoCorrection=3, you must compute the spatial gradient (e.g., from the Analysis menu) or load a previously computed spatial gradient (e.g., from the File/Parameter Map menu). You can also write your own version of computeSpatialGradient to compute whatever you want to use for inhomogeneity correction; it must save the result in the proper parameter map format in a file called "spatialGrad.mat" (see computeSpatialGradient.m for details).
Analysis steps for retinotopy data. These are the steps for analyzing retinotopy data when they were collected using both clockwise and counterclockwise wedges and both expanding and contracting rings. If only clockwise and only expanding stimuli were used, then you can skip the first 2 steps.

Time shift the tSeries (from the Analysis menu) by ~ 4-6 secs. Note that mrLoadRet shifts (to advance the phase of the response) by an integer number of frames. You will need to know the frame period (numInterleaves * TR) in order to know how many frames correspond to 4-6 secs. This time shift approximately compensates for the hemodynamic delay so that clockwise and counterclockwise data can be averaged together most effectively. It is not critical for this to perfectly compensate for the hemodynamic delay.
Flip (time-reverse) the tSeries (from the Analysis menu) for counterclockwise rotating wedge and contracting ring scans.
Average the time series (from the Analysis menu) across all repeats of the rotating wegdges, and separately across all repeats of the expanding rings.
Switch to the 'Averages' data type.
Compute the corAnal for both scans (i.e., average rotating wedges, average expanding rings)
Transform to the Gray and then onto the Flat (you'll have to open Gray and Flat windows along the way).
View the data in the Flat window, setting a reasonably high correlation threshold.
Inspect the data for reversals in the rotating wedge phase map that run orthogonal to the expanding rings phase map, and draw ROIs by hand using the ROI polygon options from the ROI menu.
Transform the resulting ROIs back to the Gray view and save them (from the File menu).
We will soon have tools for defining the ROIs on the flat map automatically given the average rotating wedge and expanding ring phase and correlation maps [stay tuned].
Analysis of blocked-alternation data sets. We have some standardized code in place for analyzing data from blocked-alternation experiments. The relevant functions are:

analyzeSessions: Loops through sessions, computing vector means for each ROI. For each ROI and each session, saves analysis structure in file under Analysis subdirectory.
analyzeConditions: Loops through conditions to compute test scan responses. Saves results in file in dir you're in when you call this function.
computeReferencePhase: Called by analyzeConditions, averaging across scans to compute the average response phase, which is needed to compute the projected amplitudes.
The steps in the analysis are as follows:
Define the ROIs (as described above).
Retrict the ROIs. For example, if your experiment used stimuli that were restricted to an annulus of the visual field, then you don't want to use the entire V1 ROI that was defined based on retinotopy measurements covering a larger part of the visual field. One recommended method for restricting the ROIs is to run a session of annulus/anti-annulus scans, average them, and choose a high correlation threshold and tight phase window to identify a subregion of each ROI that corresponds retinotopically with the stimulus annulus.
Define the sessions, conditions, and ROIs data structures (see below for an example).
The sessions structure specifies information about each scanning session.
The conditions structure specifies information about each experimental condition, that typically will have been repeated across several scanning sessions.
The ROIs data structure specifies information about each of the ROIs that you want to analyze (note again that these are the restricted ROIs).
The sessions.referenceScans field is used to compute the reference phase that we use to compensate for the hemodynamic delay. The reference phase should typically be computed using data that are separate from (and hence statistically independent from) those listed in conditions, e.g., you could use the "reference scans" (that were used for restricting the ROIs) to also compute the reference phase. An OK (but not ideal) alternative is to use all the data to compute the reference phase. The potential problem with this is that the conditions with the largest amplitude will tend to dominate in the estimate of the reference phase. Then the projected amplitudes for those conditions will be biased to be even larger. At the same time, however, the final projected amplitudes depend on an accurate estimate of the hemodynamic delay (the reference phase) so we want to throw as much as possible at it. Also, the reference phase must be computed using scans that have the same block-alternation period as the test scans.
Here's a simple example of how to define the sessions and conditions data structures, and for how to call analyzeSessions and analyzeConditions. Note, however, that in general there will be more than one session and more than one referenceScan.

homedr = 'g:/matlab/reliability';
chdir(homedr);
sessions(1).path = 'e:\mri\022701a';
sessions(1).referenceScans = 1;

conditions(1).name = 'hi';
conditions(1).sessions = ones(1,17);
conditions(1).scans = [2:2:34];
conditions(2).name = 'lo';
conditions(2).sessions = ones(1,17);
conditions(2).scans = [3:2:35];

% Note V1r7 is a restricted ROI, that was computed from V1 using a
% correlation threshold (r>0.7) and phase window [0 pi/2] on the data
% from a separate session (i.e., the data used to restrict the ROI are separate
% from those used to compute the responses).
ROIs{1} = 'V1r7';
ROIs{2} = 'V2r7';

analyzeSessions(sessions,ROIs);
analyzeConditions(conditions,sessions,ROIs);

analyzeSessions saves a file for each ROI in the Analysis subdirectory. The default filename is the ROIname, but can be specified to be something else by passing optional arguments. Each file contains an analysis structure with the following fields:
session: copy of the sessions structure
ROI: copy of the ROI structure
ROIsize: number of voxels
amps: mean amplitude for each scan in session
phs: mean phase for each scan in session
analyzeConditions saves a file for each ROI in the current directory (pwd). Again the default filename is the same as the ROIname, but can be specified with optional arguments. Each file contains a responses structure with the following fields:
name: copied from conditions
sessions: copied from conditions
scans: copied from conditions
amp: mean response amp (computed via vector mean from the individual repeats).
ph: mean response phase
projectedAmp: mean projected amplitude, computed by projecting onto mean reference phase line
repeats: vector of complex amp/ph for individual repeats
projectedAmps: vector of projected amps for individual repeats
bivariateSE:
projectedAmpSD:
projectedAmpSE:
refPhase: from computeReferencePhase.m
Event-related data analysis.
[Not written yet]

Writing Scripts:
Everything that can be done from the user interface can also be done by writing scripts that call the appropriate functions. Use hiddenInplane, hiddenVolume, hiddenGray, & hiddenFlat structures for writing such scripts. These structures are like the full blown inplane (etc) structures except that they don't have a ui substructure (i.e., no window appears on the screen). Here's an example script for computing the corAnal and transforming to the gray:
inplane = initHiddenInplane;
gray = initHiddenGray;
inplane = computeCorAnal(inplane);
gray = ip2volCorAnal(inplane,gray);
mrLoadRet itself is a script. You should feel free to make a copy of it and modify it. You could have a single copy in your home directory (make sure that it is in front of the standard version in your matlab path). Or if you like, you could have a seperate version of mrLoadRet in each and every data directory. Your personalized version might automatically open multiple windows, load the inplane anatomies, select a data type, and load the corresponding corAnal data.


Converting from mrLoadRet 2.*
Change your startup.m file to have the line 'mrLoadRet3Path' instead of the 2.* path that you used to have.
Run matlab.
Change to the directory containing the session you wan to convert (using 'cd' or 'chdir' at the matlab prompt).
Make a backup copy of the mrSESSION.mat file (I usually call it mrSESSION-25.mat), just in case. I usually do this from a unix prompt but you can do it anyway you like.
Make sure that the entire directory tree is writable by you and that you have the tSeries files online (copy them from the CD if necessary).
Type 'convert2to3' at the matlab prompt. This will convert the mrSESSION structure, create the dataTYPES structure and write them both to a new mrSESSION.mat file. This will also convert the tSeries files from .dat to .mat format.
You will get a couple warnings about converting the tSeries and trashing the existing corAnal files. Click ok on those to let it proceed.
You will also get a message about moving the Tseries folder into the Original sub-folder. You need to make an Original folder and then move the TSeries folder into Original as directed. You can do this either from a unix prompt or in windows, but it will not happen automatically. From the unix prompt type:
> cd Inplane
> mkdir Original
> mv Tseries Original/Tseries
It is a good idea to check the mrSESSION and dataTYPES structures to make sure that everything got converted properly. You might also want to change some of the fields, for example, if you had run scans with a different number of frames or different number of cycles. Eventually, we will have a nice user interface for editing these data structures but you will need to do it by hand for now. See mrLoadRet hacking for details about these data structures.
Clean up the directory tree.
Delete everything from the Volume, Gray and Flat folders except for the ROIs subdirectories. Delete the Tseries subdirectories. Delete any data files (corAnal, etc.). All of this will need to be recomputed.
If you had computed averages under mrLoadRet-2.5, then there are some Tseries/Scan* subdirectories corresponding to those averages that should be deleted. These averages will need to be recomputed.
Note that this leaves the session as if it was just recon'd. The only things that remain are the tSeries files and the ROIs. If you had computed average scans, you will have to recompute them. You will have to recompute the corAnal files and any other data analysis that you might have done. You will have to re-map the data to the Gray and Flat views. But all of this will be easily (and more quickly) recomputed in mrLoadRet-3. Do not copy or try to keep any data files (other than the tSeries and the ROIs) from mrLoadRet-2.*; they are eithere in the wrong format or incompatible with the new mrSESSION and dataTYPES structures.


Version 3.0 Release Notes
There are a few major changes that have been made in mrLoadRet 3.0, enabling new features, but also forcing you to make a bunch of changes to your own code.
MrLoadRet-2.x interplates the tSeries to the size of the inplane anatomies. This is slow (sometimes painfully slow), wasteful of disk space (e.g., the corAnal matrices are 4x larger than they need to be) and unnecessarily blurs the data. In mrLoadRet-3.0, we dumped the interpolation so that the corAnal matrices are all now the size of the cropped tSeries. The data are upsampled (via simple pixel replication) only for the purpose of display, i.e., when the pseudo-color overlay is computed. You will notice that the overlays now look like they have fat pixels instead of looking blurry. This definitely changes the results slightly because the data are no longer spatially blurred. It has been tested in head-to-head comparisons with mrLoadRet-2.5. Although the results are very similar, there is definitely a small difference. The key functions that now produce slightly different results are:
getCurDataROI
ip2VolCorAnal
ip2VolParMap
You should find it easy to convince yourself that the new results are better (less blurry), but you should check to make sure that the results you get are similar.
The corAnal matrices & parameter maps are now cell arrays, each of length equal to the number of scans. A cell of the cell arraycan be empty (if the data matrix hasn't been computed yet for that scan). Otherwise, it is an array of the appropriate size for each view type, as follows:
inplane, 3d array (y,x,slice)
volume/gray, vector of length nVoxels where nVoxels=size(volume.coords,2)
flat, 3d array (y,x,hemisphere) with NaNs where there's no data
Data types are now fully integrated. This involved making changes the directory tree structure (see above for details).

The mrSESSION structure has been changed substantially (see mrLoadRet hacking for details).

Some of the features got left behind for the time being (i.e., you should come and talk to me about how to reimplement them under mrLoadRet-3.0), including:
Analysis/CorticalMag
Analysis/EventAnalysis
Analysis/CRFAnalysis (including plotCycleAmps)
phase corAnal and complex corAnal (including computeComplexCor)

 

 

Home | mrVista WorkFlow | Scanning| Post Scanning|fMRI Data Analysis |Original Lab manual |Other Resources

Stanford Vision, Imaging Science and Technology Laboratory
Department of Psychology, Jordan Hall, Building 420
Stanford University, Stanford, CA 94305-2130

Send email to : vista@white.stanford.edu
Copyright © 2003. All Rights Reserved