Resampling fmri timeseries

[Apoorva Bhandari]

Hi, 

I'm very keen to use the methods implemented in GLMestimatesingletrial of the GLMdenoise package to get denoised estimates of single trial regression weights. I note that this function does not take a specification of the design as onset times. In my data, event onsets were not aligned to the onset of a volume acquisition. I understand that one way to address this is to temporally resample the data. 

I'm trying to figure how best to achieve such temporal resampling. I'm currently trying to build off of SPM12's slice timing correction code but it would be helpful to know if such code already exists. 

Best,
Apoorva

[Kendrick Kay]

Hi Apoorva,

We actually lately have been doing some software development work to extract the GLMestimatesingletrial.m function into a standalone repository.  It is https://github.com/kendrickkay/GLMsingle/ and we have removed the GLMestimatesingletrial.m function from the GLMdenoise repo.  (If you do a git pull, you'll see that we deleted it.)

We are still in the middle of generating example scripts and tutorials for GLMsingle, but the core function itself should be good and stable.

You are correct that the function requires the stimulus design matrix to be specified in "lock" with discrete fMRI volumes.  The way I like to set things up is to have the first brain volume correspond to time = 0 s, and so the interpretation of the corresponding row in the design matrix is that the first row corresponds to what occurs at time = 0 s.  If a stimulus appears on the screen at exactly time = 0 s, then that stimulus gets coded as a "1" in the first row.

In terms of temporal resampling, it may be a little tricky to make sure that everything is right.  You are alluding to the fact that slice time correction is a critical player in this regard, which is certainly true.  Depending on your pulse sequence flavor, you might be using multiband (simultaneous multislice) which complicates things.  Regardless, you'll want to be pretty sure you know exactly when slices are acquired.

If you want to get fancy, one idea is to do the temporal upsampling/resampling in the exact same step that corrects for slice timing differences.  We talk about this a bit in https://www.nature.com/articles/s41592-020-0941-6 (see Supp Fig 1B).   You could take your original acquisition (say, TR = 2.1 s) and just create new time-series at whatever sampling rate is most convenient -- you could even go really high like an effective sampling rate of 0.5 s (or whatever).

In terms of how to perform the temporal resampling, it's really just a use of interp1.m, i think (a built in matlab function).  There is a potentially handy utility in:

  https://github.com/kendrickkay/knkutils/blob/master/timeseries/tseriesinterp.m

Exactly when and where to do the interpolation sort of depends on exactly how you have your pipeline/data set up....

Hope this helps?

Kendrick

--

Kendrick Kay (k...@umn.edu)

Assistant Professor

Center for Magnetic Resonance Research
University of Minnesota

http://cvnlab.net

--
You received this message because you are subscribed to the Google Groups "GLMdenoise" group.
To unsubscribe from this group and stop receiving emails from it, send an email to glmdenoise+...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/glmdenoise/c7eccbb7-1e18-45c9-b45c-a9e40ee3c5b4n%40googlegroups.com.

[Apoorva Bhandari]

Thanks, Kendrick! That's very helpful. I've got the resampling working - excited to try out GLMsingle. 

[Etienne ABASSI]

Hi, I am following on this discussion as I tried to resample some fMRI data using Kenrick's code but it didn't work as I wanted.

I just wanted to know at which step of the preprocessing you were doing the resampling ?

I basically want to resample my fMRI data by going from an original TR of 1.5s to a resampled TR of 0.75s (so dividing the TR by 2)

I tried to run the code at the end of my preprocessing pipeline (fmriprep+TEDANA) and just before the GLM, then in the GLM I doubled the number of slices. Is it correct ?

If not, when would you perform it ?

I thank you very much in advance if you can give me advice !

Best,

Etienne

[Kendrick Kay]

Hi Etienne - yes, you can do the resampling at the end of your pipeline; that sound be fine.  When you say you "doubled the number of slices", I assume you mean that you doubled the number of time points (when you go from 1.5s to 0.75s)?

Kendrick

To view this discussion on the web visit https://groups.google.com/d/msgid/glmdenoise/f9b2dcf4-cfea-4da6-85ae-42fd357d3904n%40googlegroups.com.

[Etienne ABASSI]

Hi Kendrick, thank you for the quick answer.

Yes this is what I mean, I doubled the number of time points in my GLM.

I note that I am using multi-echo multi-band (ME-MB), this is why I do it after the tedana merge of the multi-echos.

So if you confirm this approach is correct, I will look at it more carefully to see if I did any mistake...

I'll keep posted on the thread in case it's useful to others.

Thank you again !

Etienne

[Kendrick Kay]

Yes -- if you are finding an error, it woudl be helpful to see what you ran and what error message you are getting.

To view this discussion on the web visit https://groups.google.com/d/msgid/glmdenoise/8dd66c20-2e68-49a5-aa7f-7f582b726b75n%40googlegroups.com.

[Etienne ABASSI]

Hi Kendrick,

I finally resolved my problem ! Here is the solution, if it helps:

For information I used the built-in "niftiread" and " niftiwrite" matlab function to load and save my 4D nifti files.

My first problem, was that using the "tseriesinterp" function after "niftiread" , the wrong dimension was used for time.

Indeed, the default is dimension 1, while I had to chose the dimension 4.

My second problem was when writing the data with  niftiwrite. I had to define the correct transformation matrix.

So I loaded the header of the nifti file using "niftiinfo" prior to the interpolation, withing the header I changed the size of the temporal dimension to match the new size.

Then I save the new interpolated 4D nifti using niftiwrite and passing the header information.

It now works correctly I think. Indeed, while running the same univariate analyses, I see a lot more false positive after the interpolation. But I guess it is due to the interpolation in itself.

I still the same cluster I'm interested in before and after the interpolation, it's why I think it works.

I have a last question regarding movement regressor.

I usually get them automatically during my preprocessing. For now, I did an interpolation on it two.

But I am wondering if I should rather do the interpolation before the preprocessing, so the movement regressor would be them extracted automatically.

Altough I am not sure there would still be correct after that...

Any thought about it ?

Thank you again !

Etienne

[Kendrick Kay]

> For information I used the built-in "niftiread" and " niftiwrite" matlab function to load and save my 4D nifti files.
>
> My first problem, was that using the "tseriesinterp" function after "niftiread" , the wrong dimension was used for time.
> Indeed, the default is dimension 1, while I had to chose the dimension 4.
>
> My second problem was when writing the data with niftiwrite. I had to define the correct transformation matrix.
> So I loaded the header of the nifti file using "niftiinfo" prior to the interpolation, withing the header I changed the size of the temporal dimension to match the new size.
> Then I save the new interpolated 4D nifti using niftiwrite and passing the header information.
>
> It now works correctly I think.

Interesting, sounds good.

> Indeed, while running the same univariate analyses, I see a lot more false positive after the interpolation. But I guess it is due to the interpolation in itself.

Yes -> It could be related to how your analysis is treating the statistical interpretation. Obviously, after upsampling, there isn't any new "fresh" sampling of the noise; instead, you are really just repeating (in a sense) the existing data. So, you need an analysis procedure that knows how to properly estimate the degrees of freedom. It appears, based on what you said, that the procedure you are using is making some sort of assumption about temporal autocorrelation and degrees of freedom...

> I have a last question regarding movement regressor.
> I usually get them automatically during my preprocessing. For now, I did an interpolation on it two.
> But I am wondering if I should rather do the interpolation before the preprocessing, so the movement regressor would be them extracted automatically.

The proper way would be to first interpolate your data, and then proceed to motion estimation.

However, my guess is that if you just interpolate your existing motion regressors, the result will be very very similar.

[Etienne ABASSI]

Great, thank you for your comments.

My aim is to then use your GLMSingle toolbox and perform MVPA. So the univariate analysis was just a sanity check.

Well noted for the motion regressors, it's all clear now !

Thank you for your help overall,

Etienne