In this blog post I want to go through the analysis pipeline of layer-dependent VASO.
I will go through the all the analysis steps that need to be done to go from raw data from the scanner to final layer profiles. The entire thing will take about 30 min (10 min analysis and 20 min explaining and browsing through data).
During the entire analysis pipeline I am using the following software packages: SPM, AFNI, and LAYNII, and gnuplot (if you want fancier plotting tools)
1.) The example data
DOWNLOAD LINK: https://goo.gl/4kMEFc (Data from other studies can also be downloaded here: https://openneuro.org/datasets/ds001547)
I use an example dataset from this study. It is a 12 min finger tapping experiment covering a small FOV of the right motor cortex.
The imaging parameters are:
- 0.74×0.74×1.25 mm
- TR = 1.5 s (For BOLD and VASO each), resulting in effective resolution of 3s.
- The time series consists of an interleaved acquisition of images with and without blood nulling.
- 10 slices
- remaining parameters can be seen at this summary pdf
- The tapping task was alternating between 30 sec of rest and 30 sec or tapping using this psychopy program.
- The FOV is double tilted
2.) Motion correction
- A manual mask is defined to avoid faulty motion estimation due to variable distortions outside the FOV.
- The non-steady-state images (the first four time points) are overwritten with steady-state images.
- The SS-SI VASO sequence acquires images with and without blood nulling interleaved (see also this blog post). Thus, the motion correction is applied separately for BOLD and VASO.
- Ideally, the motion parameters should be very similar for BOLD and VASO. This needs to be checked manually every time.
If you don’t want to use SPM because of potential reluctance or license issues when using matlab, there is also an AFNI version of motion correction here (note that the results of AFNI motion correction are not as ‘crisp’ as the SPM result).
3.) BOLD correction and contrast combination
- As explained in previous blog posts (see this blog post), there is a T2*-weighting in any EPI readout. Thus, the image with blood nulling is not solely containing VASO contrast. It also has some unwanted BOLD contamination that can counteract the negative VASO contrast. This BOLD contamination needs to be corrected for. E.g. with a division of the blood-nulled image with the not-nulled BOLD image.
- Since the blood-nulled and not-nulled image are acquired interleaved, the time series are temporally upsampled (in afni) and shifted with respect to each other, such that the respective contrasts refer to the same point in time.
- The division happens in the LAYNII program
LN_BOCO. - The
-trialBOCOcommand averages all the 12 trials together. - Since VASO has an inherent T1-weighting included, the time series can also be evaluated to provide a T1-weighted MP2RAGE-like contrast as an anatomical reference in EPI space.
4.) Quality measures
- Conventional quality measured are also estimated:
- tSNR: As measure of signal stability
- MEAN can depict artifacts that disappear in the noise level for single TR imaged
- Skew and kurtosis show how Gausian the noise distribution is. When the signal is not Gausian, it could be an indication of EPI ghosting or phase inconsistency artefacts.
- Auto-correlation can also be valuable e.g. during strong head motion.
Note, this Tutorial was using the pre-release LayNii version v1.0.0. More recent LayNii versions, have a slightly updated argument structure (e.g. -input instead of -timeseries). See help function of LN_SKEW.
5.) Functional activity
- One way to estimate the signal change is to just subtract the signal during activity with the signal during rest.
- Another way is to use a GLM implemented in all the mayor software packages.
6.) Layering
The layering has been previously already described in this blog post.
-
- The first step is to upscale the data to allow smooth layers
- Then the borderlines can be drawn (e.g. in Fsleyes)
- The layers can then be calculated with
LN_GROW_LAYERS
All LAYNII program executions are tested for version v.1.5.6 and will probably work for later versions too.
layer fMRI blog 
Hi, Renzo, this is super helpful and I learn a lot. I have a simple question. Did you also do distortion correlation in some step? or maybe I miss something?
LikeLike
Hi Ruyuan,
Yes, you are right. I did not perform any distortion correction. I found that the additional reampling during the distortion correction can result in signal blurring (https://layerfmri.com/unwanted-spatial-blurring-during-resampling/).
Here, the “anatomical” reference is taken from the same functional data that are used to quantify the signal change. So the ‘anatomical’ data are in the same distorted space as the functional data. Thus, in this case there is no distortion correction necessary.
LikeLike
Hi Renzo, this example is really good. When I perform the same pipeline on the layer V1 dataset shared on the OpenNeuro, I had a warning below
If you are performing spatial transformations on an oblique dset,
such as BOLD.nii,
or viewing/combining it with volumes of differing obliquity,
you should consider running:
3dWarp -deoblique
on this and other oblique datasets in the same session.
See 3dWarp -help for details.
++ Oblique dataset:BOLD.nii is 14.000014 degrees from plumb.
Do we need an alignment? If so, what we should align to?
LikeLike
Hi Yameng,
This this “warning” is absolutely normal and can be ignored.
It is partly comming from the fact that afni wasn’t originally developed for tilted slices.
You do not need to realign the data. The corresponding spatial resampling would only result in unwanted spatial blurring.
Best, Renzo
LikeLike
I see. Thanks for your quick and helpful reply. I also sent you an email earlier today. Would you share the V1, LGN, and V5/MT masks you used in the 2019 ISMRM abstract? Thanks in advance.
LikeLike
Hi Renzo, I have a quick question about the V1 layer data shared on the OpenNeuro. Is the TR 2.5 sec for both VASO and BOLD? But I saw on the slides from ISMRM or OHBM, the TR for VASO is 2.73 sec and for BOLD is 2.07 sec. A little confused.
LikeLike
Hi,
As explained here, https://youtu.be/KwX2rscnOxY?t=307, the sequence is using an asymmetric TR. every odd imgage has a TR of 2.73 and every even emage had a TR of 2.07, plus 200ms dead time at the end of every pair.
Since BIDS does capture this flexibility, OpenNeuro is $#%$ strict about BIDS, we needed to settle on one single value. So the TR in the header is the average TR.
LikeLike
Thanks for the reply. As you explained, in the preprocessing pipeline, we should set 2.5 sec for VASO and BOLD, right?
LikeLike
Hi Renzo,
Thanks for your tutorial.
In the parameter protocol, TR=1648.80ms. Why is TR=1.5 s in data processing? I am looking forward to your reply. Thank you!
LikeLike
Hi Wang,
Yes, this is confusing
I believe it can have multiple reasons.
1.) The TR in the Dicom header cannot deal with the alternating TRs of interleaved BOLD and VASO acquisitions.
2.) There are small missmatched of the desired TR and the actual TR
3.) It can be more straightforward to calculate the task blocks in seconds
See more explanations on the TR weirdness in VASO sequences here: https://layerfmri.com/2017/11/26/ss-si-vaso-sequence-manual/#TR_bookmark
LikeLike
Hi Renzo,
Thanks for the excellent tutorial!
I’ve went through the pipeline you’ve suggested in the video, but I’ve found that my image looks bit different from yours(Especially, scale of upsampled BOLD image and T1-weighted image). I found that your data before motion correction(Smagn.nii) has no temporal linear drift which is quite different from mine. So, I am wondering if you’ve went though some additional processing steps in during making the Smagn.nii file.
I am really appreciating your fascinating works! Thank you!
LikeLike
Hi Jungwoo,
Thanks for checking out the analysis pipeline.
It is indeed a bit strange that the you do not see the same drifts. Are you referring to the same data set that you downloaded from the link?
The Smagn file represents the raw data as they come from the scanner. They are reconstructed from the k-space data with the vendor provided pipeliend and exported as dicoms. The only “processing step” that I did was to convert them to nii.
Usually, the only case, when I apply drift correction is during the functional GLM analysis.
In case you want to apply a drift correction on your data, I would recommend the following AFNI program:
3dDetrend -prefix detrended.nii -polort 2 inputfile.nii -overwrite, where 1 is the order of the polynomial (linear only).
Very best regards,
Renzo
LikeLike