Model-Free FMRI Analysis Practical

Practical Overview

MELODIC - Single Session Analysis 

cd ~/fsl_course_data/melodic

Melodic is the tool in FSL that decomposes single FMRI data sets into time-courses and spatial maps. For group FMRI data analysis MELODIC also estimates vectors which describe the relative strength of activation in each of the single FMRI data sets. There is a simple GUI (type Melodic or Melodic_gui at a command line) which has been structured to be similar to the FEAT GUI:

The GUI will allow you to change some of the most common default options. To get full control, use the command line version: 

melodic --help

Note however, that the melodic command-line program will not do all of the preprocessing for you that the GUI-based scripting will do.

As a first example we will re-analyse the audio-visual data set that we used in the FEAT practical.

cd ~/fsl_course_data/melodic/av

and bring up the Melodic GUI (type Melodic or Melodic_gui):

Press Select 4D data and choose the av.nii.gz data set as input. This data is a cut-down dataset (only a few slices) and will bring up a warning relating to the small FOV - just ignore this warning by pressing OK. The default TR setting of 3.0 is already correct for this dataset and does not need to be changed, but in general it needs to be changed to the correct value.

Change the high-pass filter to 90s.

Select the Pre-Stats tab.

Low-frequency drifts and motion in the data can adversely affect the decomposition. In most cases, you would want to motion-correct the data and remove these drifts first. This can be done from within the Melodic GUI Pre-stats section.

The data set av.nii.gz has not been corrected for head motion; so we will leave this setting switched on.

In order to speed up calculation the data only contains 5 axial slices - in cases like this you should switch off Bet.

Spatial smoothing is not in general required for ICA, though it helps in reducing the number of spurious components; you can lower the FWHM to a small value (of the order of the voxel dimensions, 3mm in this case).

Select the Registration tab. By default MELODIC will register the middle timepoint image from the 4D FMRI input data directly to the standard space template (this is saved as example_func in the output directory). We recommend in general turning on the Main structural image option so that the lowres FMRI image is first registered to a highres structural image from the same subject. For now we will leave the non-linear registration off in order to save time, however this is generally recommended and works like it does in FEAT (requiring both brain-extracted and original structural images, whereas here we just need the brain-extracted one). For this linear registration, the brain-extracted highres image is then registered to the standard space template, and then the two registrations are combined to give an example_func2standard.mat transform which will be used later to resample the FMRI stats directly into standard space.

If the 4D data here only contains a few slices, then even before registration to the highres image, it is a good idea to register example_func to a whole-head EPI image which contains basically the same slices as the 4D data.

Select Initial structural image, set the file to epiwholehead.nii.gz and set the DOF to 3 (as we only have a few slices in the timeseries data which correspond to certain slices in the epiwholehead image, we do not want to allow the image to be rotated). Set the Main structural image file to structural_brain.nii.gz with 7 DOF (note that we have already run BET on this, and have reduced the resolution of the structural image to save time on the highres to standard registration). Leave Standard image turned on with MNI152_T1_2mm_brain.nii.gz selected and set the DOF to 12. Leave the "Resampling resolution" at the default value.


The Stats / Post-Stats sections finally let you control some of the options for the decomposition.

The default setting will most probably already be set to what you would want most of the time:

By default, MELODIC will automatically estimate the number of components from the data - or you could switch this option off and then specify the number of components explicitly.

MELODIC will also by default carry out inference on the estimated maps using the mixture model approach. A threshold level of 0.5 in the case of alternative hypothesis testing means that a voxel 'survives' thresholding as soon as the probability of being in the 'active' class (as modelled by the Gamma densities) exceeds the probability of being in the 'background' noise class. This threshold level assumes that you are assigning equal 'cost' to false-positives and false-negatives. If (for example) you consider false-positives to be twice as undesirable as false-negatives then you should change this value to 0.66...

Melodic can perform post-hoc GLM analyses of the estimated ICA time courses against a GLM design matrix and also calculate contrasts. Note that here the GLM is run after ICA has identified the time courses, i.e. as a post-processing reporting step. In order to get Melodic to calculate this you will need to specify the file names for the GLM model (XX.mat) and the contrast matrix (XX.con). The data in this example has already been analysed previously in the FEAT 1 practical. Locate the relevant files (design.mat and design.con) in the .feat directory ~/fsl_course_data/fmri/av/fmri.feat and add them to the setup.

Now simply press  Go

After a few minutes Melodic will have generated the results and your terminal window will tell you where to find the web report - though, as with FEAT, MELODIC will automatically have opened a web-page report which tells you about the progress. To see the ICA report pages, click on the ICA link.

Once Melodic has finished estimating the components, the MELODIC report page will be updated to include information text about the analysis. Also, the web report will contain information about the initial PCA dimensionality reduction and estimation of the model order. Each component then has a separate HTML page that can be accessed from the list of component numbers at the top of the main report page.

Each component page shows one spatial map thresholded and rendered on top of an example functional scan followed by the relevant time-course of the ICA decomposition and the power-spectrum of the time-course. When a design matrix and contrasts have been specified in the GUI, (as we did with the design from the FEAT 1 practical) the time-series plot will also include a full model fit and a GLM regression fit table will be reported as part of the web page. If you click on the thresholded map, you can inspect the raw IC output, the probability map (derived from the Gaussian/Gamma mixture model) and the Mixture Model fit. Use "<" or ">" to view another component map.

Amongst the set of components you should find 2 spatial maps that identify visual and auditory activation patterns:

Visual




Auditory


Dimensionality reduction

In the previous example we generated ICA results from artificial data. In that case, we used the Melodic automatic dimensionality estimation to infer the number of components to extract. In this example, we will look at the effect of dimensionality reduction.

In ~/fsl_course_data/melodic/artdata you should find three Melodic folders dim15, dim33 and dim191. These contain the Melodic report output after running on the artificial data set of 196 time points when explicitly setting the dimensionality (number of extracted components) to 15, 33 (close to the number which melodic would estimate automatically) and 191 (the number of components that you would use if you wanted to keep 99.9% of the original variance). You don't need to run Melodic again!

Instead, inspect the Melodic reports for all three runs inside the webbrowser. Open a new browser window and go to

~/fsl_course_data/melodic/artdata/dim33.ica/report/00index.html

The main results are summarised below:


Here are the three ICA maps found when using automatic dimensionality estimation (33 components) overlaid onto an example epi and thresholded using the Gaussian/Gamma mixture model...




... and the un-thresholded ICA 'Z stat' maps (click on the thresholded maps to get these)





Now we show the ICA results obtained after reducing the dimensionality to 15. Out of those 15 estimated maps, only one appears to have detected any signal. Can you explain why this happened?



Next, ICA results have been calculated after reducing the dimensionality to 191 (hardly any reduction from the original 196 timepoints). For one of these maps signal detection is significantly less accurate. Can you explain why this happens?





Resting-State data analysis

During this tutorial you will be running Melodic on resting-state FMRI data. 
cd ~/fsl_course_data/melodic/resting

In this directory there is a short resting FMRI data set (100 volumes)

First, bring up the Melodic GUI and set up the GUI to run on short.nii.gz. All default options are fine, set Threshold IC maps to 0.99 and select Output full stats folder in the Post-Stats tab.

Press Go

The analysis will take ~10 minutes. When Melodic is finished, go through the web report page and identify e.g the somatosensory RSN.

Often it helps to view the ICA maps in standard space. To do this, go into the ICA directory 

cd ~/fsl_course_data/melodic/resting/short.ica/filtered_func_data.ica

and use the already calculated registration parameters to re-slice melodic_IC: 

flirt -in melodic_IC -ref $FSLDIR/data/standard/MNI152_T1_2mm -out melodic_IC_hr -applyxfm -init ../reg/example_func2standard.mat

Now, use fslview to view the transformed maps on top of the standard brain and use the fslview atlas tool to check the locality of resting activation. 

 fslview $FSLDIR/data/standard/MNI152_T1_2mm  melodic_IC_hr.nii.gz -l "Red-Yellow" -b 5,10

Data denoising

This tutorial will give you an example of how Melodic results can be used to understand problems you might encounter when analysing your data. 
cd ~/fsl_course_data/melodic/denoise

In this directory, you will find a Feat folder denoise.feat and the output from Melodic denoise.ica

Both results were generated from data obtained on a 1.5T system under a simple motor task.

First, bring up the FEAT report page in a new browser window. You should see whole slice activation in the most inferior slice - the expected motor activation does not show up strongly.

If you inspect the Melodic report (open
~/fsl_course_data/melodic/denoise/denoise.ica/report/00index.html
in a second browser window) , you will find image artefacts that show up in each slice separately. An example spatial map looks like this:

In this particular case it turned out that there was a problem with the imaging hardware.

Go through the report pages and note down the numbers of all the affected components - or just some of the components if you are running short of time.

You can now try to filter the data. For this, you will need (i) the name of the original data, (ii) the mixing matrix that defines the decomposition and (iii) a list of component numbers to remove:

In a terminal window type 

fsl_regfilt -i denoise -d ./denoise.ica/melodic_mix -o denoise_ICAfiltered -f "a,b,c,d,e,f,..."

where a,b,c etc. are the component numbers noted down earlier.
Note: You need those doublequotes so that the entire list of numbers is passed to fsl_regfilt as the argument of the -f option!

fsl_regfilt will create a new 4D image  denoise_ICAfiltered.nii.gz
This is the original data after removal of the unwanted artefactual components.

Re-run Feat using the filtered data:

This run will take a few minutes to complete. After FEAT finishes, inspect the new results... You should now be able so see an improved sensori-motor activation pattern.

Multi-Subject analysis using Tensor-ICA

This tutorial will guide you through an entire multi-subject analysis using Melodic. 
cd ~/fsl_course_data/melodic/tensor

There are 4 directories, each containing the pre-processed data from a different subject during a simple finger tapping experiment.

First, bring up the GUI:   Melodic    or    Melodic_gui     and set up all the relevant information:

Press Go. The entire run will take about 20 minutes to complete; in the meantime you can start looking at pre-calculated results of the same analysis.

Open motor.gica/report.html and have a look at the registrations. For group ICA to work, all the data needs to be registered into a common space prior to the analysis. If registration fails it is likely that the Tensor-ICA results contain a lot of components related to this mis-registration and fewer components that relate to effects of interest.

In the MELODIC output the main page shows some summary information about the decomposition: the TICA Subject/Session modes plot contains a series of boxplots, one per estimated component. The boxplots themselves show the distribution of the relative effect sizes per subject for a given component. The components are ordered in decreasing order of median effect size. The blue box shows the inter-quartile range, i.e. 50% of all subects/sessions have an effect size within this range. These boxes therefore show how consistent an effect is across the population. Effects which are very consistent and strong across the population have a high median activation (the centre red line) and a small amount of variation around that median level (i.e. a small blue box). In contrast, effects which are outlier effects (i.e. only really exist in a small number of subjects) will have most subject-mode values close to zero and a few outliers in the boxplot, shown as single red dots outside the inter-quartile range. Note that we here only have 4 data points for each boxplot, so quantities like the inter-quartile-range will be estimated quite poorly and outliers unlikely to be identified.

The second plot (PCA estimates) shows the Eigenspectrum of the data - this has been used to estimate the number of components.

Now, start looking at the third component map. Each component shows three main things: spatial map, time course and subject/session effect sizes. Both the time course (with power spectrum) and subject/session effect size plots also have associated GLM results if a design matrix and contrasts are specified when running Melodic.

The components are ranked in order of the median effect size across the population, however the best fitting component to the design is, in this case, component number three (if you want, check the F-statistic for the full model fit in the first few components). Its spatial map shows clear motor activation with a reasonably symmetrical pattern.

During the experiment subjects performed either index-finger tapping, sequential finger tapping or random finger tapping with their right hand. Here is Tdesign.png, the single-session FEAT design matrix:

In Melodic we have used the design matrix and the associated contrast matrix to test the time courses for significant correlation with the task. The time-series plot also contains the total model fit of the GLM design against the estimated time course. Underneath the powerspectrum, Melodic has generated a table of the relevant statistics which contains the GLM parameter estimates (note how the amplitude changes when under different types of task), a simple F-test of the total model fit and t-test for each of the lower-level contrasts. The time course is significantly 'active' during all three tasks. Also, the GLM fit suggests that the amplitude change during random finger tapping is larger than during sequential finger tapping (COPE 4) which itself appears to be significantly larger than the amplitude change during index-finger tapping only (COPE 5)

Beneath this GLM table the web page shows the relative signal strength for each one of the 4 inputs (subjects). The second design matrix Sdesign.mat was used to test for significant mean group activation in the same way that this is done in a higher-level FEAT analysis, see Sdesign.png:

ICA versus PCA (Optional)

In  ~/fsl_course_data/melodic/artdata/ you will find an artificial FMRI data set. Go to the directory and run Melodic from the command line: 
melodic -i artdata -v --nobet --report --Opca --Owhite

Bet has been switched off here because the data only contains 3 slices. The -v option simply turns on verbose reporting of what is going on inside melodic and the   --Opca, --Owhite  requests additional output information which we are going to use later.

With the --Opca --Owhite we requested that melodic saves the output from the intermediate PCA step. In order to compare the ICA results with results from a standard PCA analysis we now need to perform mixture-modelling on these PCA maps and generate a second report. We will therefore call melodic to run mixture modelling only: type 

melodic -i artdata --ICs=artdata.ica/melodic_pca \
  --mix=artdata.ica/melodic_dewhite -v --report -o pcaresults

Whilst this is running, read on...

The original melodic command above will have created a folder called artdata.ica which contains all the information that Melodic outputs.

The  --report   option generated a subdirectory report  inside  artdata.ica in which melodic generates web pages that show the IC maps overlayed on top of an example EPI together with the associated time-courses etc. Point your browser window at

~/fsl_course_data/melodic/artdata/artdata.ica/report/00index.html

Melodic performed mixture modelling on the IC maps - if you click on any of those maps you can inspect the details of the mixture model fit.

You should find 3 spatial maps clearly corresponding to the simulated added "activations".

Once the mixture modelling on the PCA output finishes, look at its equivalent output; there will now be a second output directory called  pcaresults  - this one contains the report directory for the PCA output. Inspect these results in a webbrowser window and try to find the three 'activation' maps, by looking through all the non-thresholded maps in the web report.

The three letters ("activations") are contained primarily in three different spatial maps in the ICA output.
Can you explain why the PCA results do not identify the same sources as the ICA analysis?