Removing nuisance PCA components from the fMRI data

Question

So in attempting to replicate analysis pipeline from Tambini & Davachi, PNAS 2013, Persistence of hippocampal multivoxel patterns into postencoding rest is related to memory I'm hoping to use PCA as a noise reduction technique in MATLAB but I have never done this before and have some technical questions. I'm coding this myself and am very unsure I'm doing it correctly, which worries me.

Basically my input are individual NII files (128x128x33x240, here first three dimensions encode voxel coordinate in space and the last one encodes time) of masked CSF/WM (Cerebrospinal fluid & white matter) activations which I'm hoping to run PCA on - these are areas of no interest in that activations here do not represent 'real' neural data and so PCA is done on them so the most robust components can be removed via regression. Then I'd like to retain the top however-many components (e.g. top 8), create an output file consisting of just those components and use that output file in a regression via the AFNI program 3dTproject. Most of my questions are technical but they are numerous so I apologize.

1. Given that the input data sets are already masked - i.e. everything outside of these areas is set to zero so that in my original 128x128x33x240 matrix, I have mostly zeros and then values at the 3 dimensions in space (at all 240 time points) where the WM/CSF resides area equal to 'activation' values - in that everything outside of the WM/CSF (white matter & cerebrospinal fluid) is equal to zero, should I mask again once in MATLAB (with logical arrays or the like)? In other words, if I input the original 128x128x33x240 image with all of the zeros will all the zeroes affect the PCA in any way?

2. Okay this is a big one for me. I have done a lot of reading on this and think I have the answer but want to confirm because of how important it is. So if I ultimately use the MATLAB pca function to run the PCA I need a 2D input. This input should be VOXELSxTIME and not TIMExVOXELS, correct? So e.g. will be 25654x240 and not 240x25654? I understand that the number of PCs will truncate at 240 regardless of the input matrix dimensions but obviously this needs to be correct in order to get the right output.

3. Okay so my current command for actually running the PCA is:

    [COEFF, SCORE, LATENT] = pca(image);
   

   COEFF outputs a 240x240 matrix. SCORE is VOXELx240. LATENT is 240x1. I know LATENT represents the eigenvalues which I will be using to choose which components to keep. COEFF is the eigenvectors? SCORE is said (in the help) to be the representation of the input image into the PC space.

   After figuring out how many PCs to keep I want to create my output file consisting of ONLY the top 8 (or whatever) PCs I'm retaining and I'm doing this by multiplying the original image (e.g. 25654x240) by the eigenvector matrix of the retained components (e.g. 240x8) so the result will be a 25654x8 matrix. This is my problem - I don't know what to do from here. I want to be able to use the PCs as a regressor and to have my output be a 3DxTime output. The exact line in the manuscript here is:

   After selecting the number of PCs, the associated temporal projection or score of each selected PC was then removed from the BOLD data in a voxelwise fashion by linear regression.

   The exact section I'm looking at is in the Supplementary Materials page 2/11 and is the bottom left paragraph (under heading, "Functional MRI Data Acquisition and Preprocessing" all the way to the end of the section).

   If my understanding is correct, each column of my matrix represents the voxel-wide representation of the 1st, 2nd, etc. PC. So what does "temporal projection" mean in this case? And how do I use this in a regression?

   I ask this because the regression programs in fMRI analysis softwares like AFNI's 3dTproject take input files of size TIMEx1 (so 240x1 in my case) which represents temporal PCA as opposed to a spatial map, which I have here. I'm assuming I have to reshape the matrix to the original size (e.g. 128x128x33x8) with the fourth dimension now being 8 instead of 240 so 128x128x33x1 is PC image 1, 128x128x33x2 is PC image 2, etc. Do I save all the 4D matrices separately? If I then use 3dTproject to regress these out, what option do I use as it is literally a single image and not a timeseries (240x1) file? Or does the 'temporal project' statement indicate that I DO need a time component to these inputs and if so how do I get that to be part of it here?

   Is singular value decomposition something I should be looking into?

This is a big question and I would be so grateful for any help. I have been trying to work on this for awhile but I feel like my understanding is not where it needs to be and I definitely don't want to mess this up.

Answer 1

Let us pause for a second and think about the purpose of this analysis. The recorded data consist of 240 time points for each of the 540672 (128x128x33) voxels. Some of this activity is presumably not relevant for the experiment (perhaps global modulations over the whole brain or correlated measurment noise), and the idea is:

1. to identify this "nuisance" activity by looking at a subset of voxels corresponding to white matter or cerebrospinal fluid, i.e. not to neurons,
2. and then to remove this nuisance activity from all voxels.

So what you want to do is to look at the nuisance ROIs (let's say there is 100000 of them) and extract some leading PCs (let's say 8). This means that voxels are your dimensions/features and time points are your data points. You want to do PCA and to reduce the dimensionality from 100000 to 8.

This allows to answer your questions:

• Question 2. The input to pca() should be of 240x540672 size because data points are supposed to be in rows and features in columns.

• Question 1. If the non-nuisance voxels are all zeroed out in the matrix, then you don't need to remove them. These will be variables with zero variance and they will not contribute anything to the PCA results. (But removing them can accelerate the computation quite a bit.)

• Question 3. If nuisanceData is an array of 240x540672 size and you run

    [COEFF, SCORE, LATENT] = pca(nuisanceData);
  

  You will get principal components in SCORE. You only need first 8 columns of it: SCORE(:,1:8). This is the matrix of 240x8 size and it contains nuisance activity patterns. Now for each voxel in the full data you need to regress its activity onto these 8 columns and only keep the residuals. If properly coded, this can be one line of Matlab code.

  When the paper talks about "temporal projection or score of each selected PC" they mean one column of SCORE matrix. It is a time series (consists of 240 time points), hence probably why they call it "temporal projection", but I agree that this term is somewhat confusing.

As a last remark: running pca() on 240x540672 matrix might take quite long. You can perform PCA much faster by coding it yourself via eigen-decomposition of the Gram matrix.