Beauchamp:RestingState

From OpenWetWare
Jump to navigationJump to search
Brain picture
Beauchamp Lab Notebook






How to acquire a resting state scan

The idea of a resting state experiment is pretty straightforward-- you put the participant in the magnet and don't ask them to do much. There are variations on "don't do much," including "Keep your eyes open and stay awake," "Please watch the crosshair/sticker on the inside of the bore and don't think about anything in particular," and "ok, we're all finished with the task now." In general, you want your participant to be awake if possible, and to the greatest extent possible, not moving.

For the acquisition itself, one can generally copy their fMRI task acquisition. (In general) You will use an epi scan with TR in the range of 1-3 seconds (with shorter generally being better), a TE in the range of 25-40 seconds, and voxel size 2-4 mm (with smaller being better). It's common to use accelerations (multi-band/simultaneous multi-slice and GRAPPA) to accelerate the sequence so you can get shorter TRs and higher resolution. Some common approaches are GRAPPA=2 and MB=3 or GRAPPA=0 and MB=6 or 8. Longer acquisitions are better. You should be aiming for at least 6 minutes for healthy adults who are going to stay still (8 minutes for wigglier individuals), and 10-15 minutes is much better. Some groups acquire much longer (30-120 minutes, sometimes over multiple sessions) and this gets some really cool stuff, but that's a bit of a commitment. Two 6- to 8-minute runs is a good place to start. Breaking the acquisition into two runs is useful for a few reasons. First, you can remind your participant to stay still and stay awake at the half-way point. You can also check over the first half the data here and make sure they're staying still enough (and possibly make the choice to re-acquire if it's not up to par). Also, you can acquire the two runs with opposite phase encode directions (that is, "blip up/blip down," which is commonly labeled AP/PA) and then the distortion differences between the two can be used to correct for susceptibility artifacts generated by the sinuses. That is useful for later registration to your anatomical images (which have less distortion).

Overall approach

The basic idea of a resting state analysis is to identify areas of the brain that co-vary over time. Areas that consistently co-vary together are considered part of the same "network." Most of these networks are well-characterized and have names. (Here's a good reference, see figures 11 and 12). There are a few ways to identify these in your own data, the most common being seed-based approaches or ICA. Let's start by considering a seed-based approach because it's simpler to think about. You'll pick your "seed," which is known to be part of a network. You find the timecourse within that seed and correlate it to everything else in the brain, assigning a correlation value to each voxel. Now, you have a map of correlation values. You can take the average of these (over your subjects, see Creating Volume Average Datasets with AFNI) to find out what's in the same network as the seed. You can also compare the correlations across subjects or conditions (using GLMs, ANOVAs, etc). NB: It's common to find that the parts of the map that change across participants/sessions are not "in the network," because the network is the set of regions that are consistently highly correlated, and thus unlikely to vary across people.

IMPORTANT: Noise regression is pretty much the most important part of resting state analyses. When you calculate your map of correlations, there is noise (often, shared noise) in both the seed region and in the places you are correlating the seed to. Shared noise inflates correlation values, and noise that varies across the brain biases them low. One MUST be careful about how they are considering these sources of noise. Especially when you are comparing across sessions/people and not across brain regions. (When you're comparing across brain regions, which is to say, determining if region A or region B is more correlated to the seed), the shared noise is probably about the same for both regions, and it should mostly come out in the wash when you average over people. If you're comparing across people, you're asking if region A in person 1 is more or less correlated than region A in person 2, and each of these has different sources of noise that you won't be able to average away. This is, of course, especially problematic if all of Group I and Group II have differences in noise.

Example afni_proc.py script - some reasonable choices for noise regression

As with so much in MRI, there is no right way to process your data (although there are some wrong ways). Here's an example of what I've used before.

Start with running FreeSurfer and convert to SUMA format. $sid here is your subject identifier. 

 recon-all -all -s $sid -i $path_to_T1_volume
 @SUMA_Make_Spec_FS -NIFTI -sid $sid

Then, make masks of the ventricles and white matter

(In the SUMA folder-- these numbers come from the labels that are output by FreeSurfer) 

 3dcalc -a aparc+aseg.nii -datum byte -prefix vent.nii -expr 'amongst(a,4,43)' 
 3dcalc -a aparc+aseg.nii -datum byte -prefix white.nii -expr 'amongst(a,2,7,41,46,251,252,253,254,255)'

Now, set up a noise regression afni_proc.py script. In this example, we're

  1. volume registering the resting state dataset (that is, aligning each volume to all the rest). During this step, it's also noting how much each volume moved and keeping both the mean motion and the derivative in each direction. We'll use these motion traces to regress out a bit later. The idea here is that if a subject moves a small amount, and a given voxel is on the edge of something brighter and something darker, then regressing out movement in the direction perpendicular to that border will compensate for this partial volume averaging. If the movement is too big, this won't hold (assuming it ever did), but that's why we censor for motion as well. We are keeping the derivative, too, because (especially with interleaved acquisitions, as most fMRI acquisitions are), there are effects that persist past one TR.
  2. aligning the resting state dataset to the T1 using the LPC+ZZ cost function, which will then be registered to the TLRC atlas (I also start by aligning their centers of mass, but if they were acquired in the same session, that's probably not necessary, but it also shouldn't hurt)
  3. bringing the two anatomical masks (ventricles and white matter) I just made along for that ride
  4. eroding the ventricles and white matter (that is, removing the voxels that are on the border so that what's left is definitely inside the mask), and calculating the first PCA component.
  5. blurring by 5mm. In this example, the voxels are 3.5mm. My rule of thumb is to blur by 1.5 voxels.
  6. Using anaticor, which takes a very blurred version of the white matter and creates a per-voxel regressor. The idea here is that BOLD signal is in the grey matter, but noise will be shared between grey matter and nearby white matter. This is finding the signal from the nearby white matter to remove it.
  7. censoring: (or, scrubbing) If there was a volume where either the participant moved more than .2mm or 10% of the voxels are temporal outliers. 0.2 mm is a fairly arbitrary cutoff. I would set this lower for smaller voxels and shorter TRs, but in the end, it's really a trade off for how many time points you will lose (and maybe how many participants you'll have to throw out). It's best to check this decision against a histogram of how *your* subjects did, and make sure you're not losing more data than you can afford, introducing bias by selecting against too many subjects or certain types of subjects, or losing power such that your findings are unstable. For a general point of reference, in most studies I've processed, I've lost ~10% of subjects to motion. YMMV. If you're comparing between two groups, make sure to check after censoring/excluding subjects that your two groups have comparable leftover motion and degrees of freedom after censoring.
  8. regressing: 3 Legendre polynomials, the first PCA component of the ventricles and white matter, the motion trace and derivative of that, the ANATICOR voxel-wise regressor.
 afni_proc.py -subj_id ${n} \
  -blocks align tlrc volreg mask blur regress \
  -copy_anat ${T1_path}/MECIG_${n}_SurfVol.nii \
  -anat_follower_ROI vent epi $T1_path/vent.nii               \
  -anat_follower_ROI wm epi $T1_path/wm.nii \
  -anat_follower_erode vent wm \
  -dsets RS_${n}.nii -align_opts_aea -cost lpc+ZZ \
  -volreg_align_e2a  -volreg_tlrc_warp  \
  -tlrc_opts_at -init_xform CENTER \
  -blur_size 5 -regress_polort 3 -regress_censor_outliers .1         \
  -regress_censor_motion .2 -regress_apply_mot_types demean deriv      \
  -regress_ROI_PC vent 1 -regress_ROI_PC wm 1 -regress_anaticor_fast   \
  -regress_anaticor_label wm

My output of interest from this is the errts.fanaticor dataset. This is a "cleaned up" version of the whole brain time course. This is, overall, a pretty conservative processing stream, meaning that we're regressing out a lot. Remember that each regressor is only an estimate of noise, and it will be correlated to some degree with real neural fluctuations (which you will then remove from your signal). Consider your degrees of freedom... don't regress and censor away all of your statistical power. If you're running a between-subjects contrast, it's also not a bad idea to check to see if any of the statistical maps to these regressors differ between groups. (So, for example that I'm completely making up, lets say you have one group with white matter disease, and the pattern of blood flow in their white matter is very inconsistent. When you take the first pca component in the white matter, you capture less variance than in your healthy controls, and have thus introduced bias by this sort of "cleaning up" of your data.)

Some reasonable options not included in this particular example:

  1. If you want to do a CompCor approach, you could use something like 5 pca components instead of just 1. But again, each one will be associated with some amount of "real" (that is, neural) variance, and if you take too many components, you might just start to regress out your resting state networks with the largest real estate. Also, since you're fitting each one separately, you're going to start to look more and more like a global signal. And if you're going to do that, there are much more computationally effective ways to tick off your reviewers. ;-)
  2. It may not be so necessary to regress out both the PCA of the white matter and ANATICOR. In most cases and brain regions, they're probably very similar. And if they're not, you should probably know why...
  3. You could register and blur on the surface rather than in the volume. To change that, drop the "tlrc" and "blur" steps. You'll end up with an un-blurred errts time course, that you can sample to your surface and then blur along the cortical ribbon. There are surface examples in the afni_proc.py help file, but they don't seem to work when you're using ANATICOR. Also, consider-- when you process on the surface, you've effectively applied a grey-matter mask. This is fine and dandy for interpreting your results, but it makes it very difficult to check for funny business happening during your processing (like massive correlations outside the brain, and the like). So, it's probably good practice to run a volume analysis anyway, just to check things over.
  4. Bandpass filtering. A common step is to exclude all frequencies lower (that is, slower-moving) than 0.01 Hz (they're indistinguishable from scanner drift) and any frequency above either 0.08 or 0.1 Hz. The rationale here is that the hemodynamic response is only so fast, and anything happening faster than that is just noise. You will probably have more correlated data (that is, higher connectivity) if you do this step, but you'll also lose a lot of degrees of freedom. If your TR was 2s, than the highest frequency in your data is (1/(2*TR)), or 0.25 Hz. Throwing away all data above 0.10 Hz means getting rid of 60% of your data. (!!!!) This is even worse if you have a shorter TR. (If your TR is 800ms, you are sampling up to .625Hz, and you throw out 84% of your data this way.) If you want to do it, you should use 3dTfitter or similar to bandpass at the same time as you do the rest of the regressions. (So, really, instead of bandpassing in the frequency domain, you're regressing out a series of sines and cosines). What you definitely definitely do not want to do is (1) regress with censoring, then bandpass (because you've now broken the time axis), or (2) bandpass, then regress non-filtered regressors (because now you're re-introducing high-frequencies, but only if you got them from the noise estimates). So, simultaneously handling these is best.

You may have noticed that there's no output correlation map here. That's because (1) I want to do some QA first, and (2) from the end of this step, we can do either seed analysis or ICA. So I'll come back to add the seed after we check that this much looks reasonable.

Quick QA

At this point, I use InstaCorr to check my data. (OK, there's a lot more QA to be done, but this is the resting-state-analysis-specific stuff. afni_proc.py will generate a QA script to walk you through a lot of the rest, and I am by no means suggesting you skip that part. It's important). If you ran the afni_proc.py script above (and the resulting script), you'll have a $sid.results folder than includes anat.final and errts.$sid.fanaticor. Load those up into afni. 

 afni anat_final.01+tlrc. errts.01.fanaticor+tlrc.

Click "Define Overlay"

And now, "SetUp ICorr"

Choose dataset. This will be the errts file in this case. You could also do this with a raw dataset (that's a great way to check for artifacts, actually). This should be a 3d+time dataset. In this example, I only told afni to read in two files, so the choice is pretty easy.

Set up the other options. We've already processed these data, so there's no need to do more. -- No need to blur. (Unless you didn't blur because you're going to sample to the surface next, then set this to be ~1.5 times your voxel size). -- I don't like to mask at this stage, because I'd like to see if there's weird stuff happening outside the brain. -- No need to bandpass -- But, we do want to have a seed that's bigger than a single voxel. So, check "Misc Opts," and bump up the radius to something larger, maybe 5. You can also turn off the polort step (we did that already in the processing.)

Use these settings, and then click "Setup+Keep" (or "Setup+Quit," if you're feeling sure of yourself.) Give it a sec to calculate things.

Now, we're going to start by finding the DMN. It's one of the most robust resting state networks, so if you can't find it in your processed data, there's a good chance something has gone horribly, horribly wrong. I find this by choosing a sagittal slice just off the midline (in about the center of that line of grey matter), and finding where the two marked sulci intersect. They have names, I don't remember what they are. But, it's a wiki page, so someone could add them. ctrl-shift-click.

Huzzah!! Now you should have a colored map over your anatomical set. And hopefully, it looks a lot like this:

There should be a big blob near where you clicked, and something else in the front. There shouldn't be that much outside the head. If you look at other views, it should be largely symmetric, and you should be able to clearly see 4 big blobs, which are generally visible in one axial slice. There should be a big blob where you clicked (there should always be a big blob where you put your seed- it should be strongly correlated to itself.) There should be two near the angular gyri, and one in the medial frontal lobes.

You can now ctrl-shift-click and drag. If you run along an axial or coronal slice, you should be able to see some largely symmetric blobs forming that more/less trace grey matter boundaries. If you pick out some areas of interest to you (motor system, visual system, etc), you should be able to see the rest of the network light up.

Seed-based approach

Now that you have clean data and you've checked it over, let's do something useful with it. One thing we can do is a seed-based approach.

In a seed-based approach, you pick one region of the brain (but see note below) and call it special. Generally, you get this region from either published MNI/TLRC coordinates and a sphere of some size around it, or from an atlas with regions marked. My preferred method is to use regions (or parts of regions) from FreeSurfer parcellations or other surface-based parcellations. (see Cortical Surface models overview and Using the HCP Atlas Area Labels from Glasser et al. in AFNI/SUMA). I prefer this because it limits the region to the grey matter. This may not be a huge deal if you're talking about young, healthy controls with nice plump cortical ribbons, but if you are doing, for example, group comparisons in older adults with cortical atrophy, a sphere on PCC coordinates may actually be mostly CSF (or even worse, more CSF in one group than the other), and that's not really what you care about (probably.) We just spent all that time regressing out noise, it would be a shame to not take similar care in finding our seed time series. I also prefer to calculate the seed time course from my cleaned up dataset. Not everyone does it that way, but that's my personal preference. The reasonable choices are:

  • Do motion correction, then calculate all the region-based time series (seed, ventricles, white matter) and put them all in one big happy regression along with your noise regressors.
  • (Do motion correction and) Do all of your data cleanup, calculate the seed time series, and then regress the "clean" seed against a "clean" volume dataset. BUT, make sure you're accounting for all those degrees of freedom you gave up earlier. So... you may still end up throwing them all in together in one big happy regression. You can, alternatively add "dummy" regressors that are orthogonal to your dataset, but still account for degrees of freedom. In practice, that looks like doing the same 3dDeconvolve command you used to get your "clean" data, but the input *is* the "clean" data, and you add the seed regressor step. That should result in a lot of beta weights of zero (hence, "dummy" regressors).


Note: Your seed doesn't actually have to be a region. Really, it's a timecourse. More often than not, this timecourse is the average within a region (the seed region). But, it could be the average from a whole network, or the average of a region convolved with a task design (that is, PPI), or a time course from a seed that has been filtered/processed somehow. All the remaining principles apply. People talk about it as a seed region, but they're really talking about an average seed timecourse, and revealing their lack of creativity by simply taking an average. ;-)


For those familiar with task-based fMRI, you can also think of this as the "task" being whatever a particular piece of the brain is doing (say, precuneus/PCC, a node in the DMN). You're computing the "activation" map where the task is "do whatever the PCC does."

Here, we'll walk through the example of using a DMN seed from the Glasser segmentation. I'm assuming here that after you ran FreeSurfer, and @Make_Spec_FS, and that you've copied the std.141 version of the Glasser segmentation over. The seed I'm going to use is labeled #34 in the Glasser segmentation, and has label "d23ab," which means Area dorsal 23 a+b. This is in the precuneus/posterior cingulate, a classic seed for finding the default mode network.

Put the Atlas into your own subject's functional space

Most of this process is described in Using the HCP Atlas Area Labels from Glasser et al. in AFNI/SUMA.

We're going to take the Glasser atlas in std.141 and put it on your own subject's space, as described on the other page. This time, since our goal is to convert to functional space with a lower resolution, I'm going to start by blurring my cortical ribbon just a bit so that more voxels are included. The ?h.ribbon.nii files are located in the SUMA folder. 

 foreach hemi (lh rh)
   3dcalc -a ${hemi}.ribbon.nii -datum float -prefix ${hemi}.fl.ribbon.nii -expr 'a'
   3dmerge -1blur_fwhm 2 -doall -datum float -prefix ${hemi}.blur.ribbon.nii ${hemi}.fl.ribbon.nii
 end

What this has done is taken the ribbon mask (which is all ones and zeros and stored as a byte), converted it to float, and blurred by 2mm. 

 foreach hemi (lh rh)
   @surf_to_vol_spackle -spec std.141.${subj}_${hemi}.spec -surfA smoothwm -surfB pial -maskset $hemi.blur.ribbon.nii’<0.25..2>’ -surfset std.141.{$hemi}.HCP.annot.niml.dset -mode -prefix HCP.volume_{$hemi}
 end

This took the blurred ribbon, and set a threshold on the blurring at 0.25. The upper limit doesn't matter that much... nothing in the dataset should be > 1. You could set this at 1 or 5000 and you should get the same answer. Now, since we want a bilateral seed, we're going to put to two hemispheres together. If you only wanted a seed on the left side, you wouldn't add the two together (because the number labels are the same and it would be a pain to tell them apart now), and you also wouldn't really have to bother computing both hemispheres, I suppose. I'm also going to specify the seed now instead of converting all of the atlas regions. 

 3dcalc -datum byte -prefix HCP.d23ab.both.nii -a HCP.volume_lh.nii'<34>' -b HCP.volume_rh.nii'<34>' -expr "max(a,b)"

You now have a dataset with all the Glasser regions that matches your subject's particular anatomy. This is lovely, but isn't quite ready yet, because we still need the voxels in the atlas to match the voxels in the functional data, and they need to be in TLRC template space, because that's what we did (in this example) to the functional data. 

 3dAllineate -base ${subj}_SurfVol_ns+tlrc -input HCP.d23ab.both.nii -1Dmatrix_apply warp.anat.Xat.1D -final NN -prefix HCP.d23ab.tlrc.nii
 3dresample -master errts.${subj}.fanaticor+tlrc -prefix HCP.d23ab.res.nii -inset HCP.d23ab.tlrc.nii

If you did the above, you can skip the next section

As an alternative to the atlas region, you could just use a sphere

If you didn't want to do all this atlas stuff, an (easier?) approach is to just use a sphere centered on the coordinate of your choosing. I mentioned above why this isn't my favorite approach, but it's still an approach you can use. To do that, find a reference somewhere that lists the coordinates of the seed you want IN THE TEMPLATE SPACE YOU ARE USING. It does no good to use MNI coordinates if all of your subjects are in TLRC space, or vice versa. Make sure you get these straight. Double check it on your anatomical.

Once you have the coordinates, generating the mask is pretty straightforward. To take an example almost exactly from the afni help, let's assume the coordinates you found were (20, 30, 70) and you want a sphere of radius 3 (sqrt(9)). Compute a sphere on the grid of your errts file: 

 3dcalc -a errts.$sid.fanaticor  -expr 'step(9-(x-20)*(x-20)-(y-30)*(y-30)-(z-70)*(z-70))' -prefix seed_mask

You'll want to view the sphere as an overlay on anat.final, just to be sure it's in the right spot, and then proceed with the steps below

Get a time course to regress and regress it

Now that we've defined our seed region in the correct space, sample the time course from it. 

 3dmaskave -q -mask HCP.d23ab.res.nii > seed.d23ab.txt

And finally, you can finish up the most exciting part of the analysis: computing a correlation value! 

 3dTcorr1D -pearson -prefix Seed_corr_d23ab.nii errts.${sid}.fanaticor+tlrc. seed.d23ab.txt

If you view this, you'll now have a map of correlation values to your seed.

Computing accurate statistics for my seed time course

You could compute the statistics of the seed to the dataset. Normally, to move forward with your analysis, you're only interested in the correlation values, not the p-value or t-stat for the correlation. So, for most people, this step isn't really necessary. But, if this is something you want, you need to account for the degrees of freedom that you lost due to the nuisance regression. Here's an example of how to do that:

I took the 3dDeconvolve command that was in the proc.${sub} script and modified it in a few ways.

  1. The input is errts.$subj.fanaticor+tlrc, rather than the motion-corrected file.
  2. Not strictly necessary, but I did it anyway-- remove the "-stim_base #" flags, so that you can confirm that all of these fits are ~0 and these are dummy variables (to preserve degrees of freedom).
  3. Add the seed timecourse as stimulus 13.
  4. Don't tell it to stop after creating the 1D file (remove -X1D_stop)
  5. add -rout so that it computes the correlation
  6. change the names of the output files. I just added ".seed"
 3dDeconvolve -input errts.${subj}.fanaticor                           \
     -censor censor_${subj}_combined_2.1D                                   \
     -ortvec ROIPC.vent.1D ROIPC.vent                                       \
     -ortvec ROIPC.wm.1D ROIPC.wm                                           \
     -polort 3                                                              \
     -num_stimts 13                                                         \
     -stim_file 1 motion_demean.1D'[0]' -stim_label 1 roll_01  \
   -stim_file 2 motion_demean.1D'[1]' -stim_label 2 pitch_01 \
   -stim_file 3 motion_demean.1D'[2]' -stim_label 3 yaw_01   \
   -stim_file 4 motion_demean.1D'[3]' -stim_label 4 dS_01    \
   -stim_file 5 motion_demean.1D'[4]' -stim_label 5 dL_01    \
   -stim_file 6 motion_demean.1D'[5]' -stim_label 6 dP_01    \
   -stim_file 7 motion_deriv.1D'[0]' -stim_label 7 roll_02   \
   -stim_file 8 motion_deriv.1D'[1]' -stim_label 8 pitch_02  \
   -stim_file 9 motion_deriv.1D'[2]' -stim_label 9 yaw_02    \
   -stim_file 10 motion_deriv.1D'[3]' -stim_label 10 dS_02  \
   -stim_file 11 motion_deriv.1D'[4]' -stim_label 11 dL_02  \
   -stim_file 12 motion_deriv.1D'[5]' -stim_label 12 dP_02  \
   -stim_file 13 seed.d23ab.txt -stim_label 13 d23ab \
   -fout -tout -rout -x1D X.xmat.1D -xjpeg X.jpg                                \
   -x1D_uncensored X.nocensor.xmat.1D                                     \
   -fitts fitts.$subj.seed                                                     \
   -errts errts.${subj}.seed                                                   \
   -bucket stats.$subj.seed

If this worked correctly, you should see a series of scary-looking warnings:

These are good. What they mean is that there are no voxels that are significantly fit by the noise regressors you've included. That's good because you already regressed them out. They're only there so that 3dDeconvolve knows how many degrees of freedom you really have. Note that you don't get the same warning for your seed regression, because that's new variation.

Now, you can view the stats.$sid.seed file the same as you would a functional dataset.

As a sanity check, at this point, your InstaCorr, stats.$sid.seed, and 3dTcorr output should all look about the same.

More than one subject

In almost all cases (nowadays) you're going to be collecting more than just one subject. You're either combining all your participants together to learn something about your seed, or you're comparing across members of your study cohort to understand something about them. If the first, you're going to follow the procedures listed here: Creating Volume Average Datasets with AFNI

If the second, we'll walk through an example using 3dMVM.

3dMVM is a powerful tool, but it's a little finicky. Spaces, quotes, and other punctuation can make the script fail, or have it run a different regression than you were looking for. And then troubleshooting can be difficult because who really notices how many spaces there are?!? Thankfully, the folks at AFNI are aware of this, and have a tool to help you build a 3dMVM script. Much like uber_subject.py to help you set up your afni_proc.py script, you can use 3dMVM_validator to help set up your 3dMVM script. All you need to start is a data table.

You can use whatever means necessary to set up the data table. For whatever it's worth (possibly very little), I find it easiest to start with an excel spreadsheet.

  1. Each subject gets a row, each data value gets a column.
  2. The first row is headers for the columns.
  3. The first column is "Subj." The second-to-last column will be "InputFile," and the very last column is the character "\."
  4. Fill in stuff like age, handedness, group, etc. etc. Any value that applies to a whole participant, fill in now.
  5. In the column for "InputFile," put the whole path to a single brick of imaging data. Don't forget that you can specify a single subbrick using the "[]" characters (e.g., MyStats.nii'[2]' to use the 3rd subbrick. That's right, third. The first one is '[0]').
  6. If you have more than one imaging "result" per subject, (for example, you have Task A and Task B),...
    1. Add a column for "Task" (or whatever you'd like to call it)
    2. Put the same value for all subjects in this new column ("TaskA")
    3. Copy your whole table (minus the headers) and paste below, now replace "TaskA" with "TaskB"
    4. Make sure you change the "InputFile" to point to the right spot.
  7. Check that there are no spaces in any of your cells
  8. Delete the "\" from the last row.
  9. Copy this into a text editor, or save as .csv.

Here's a really mundane example, which we'll pretend is saved to the file ILL_data.txt: 

Subj Age Task InputFile \
Fred 55 TaskA Fred.stats.nii'[0]' \
Ethel 53 TaskA Ethel.stats.nii'[0]' \
Lucy 38 Task A Lucy.stats.nii'[0]' \
Rick 40 Task A Rick.stats.nii'[0]' \
Fred 55 TaskB Fred.stats.nii'[3]' \
Ethel 53 TaskB Ethel.stats.nii'[3]' \
Lucy 38 Task B Lucy.stats.nii'[3]' \
Rick 40 Task B Rick.stats.nii'[3]'

Now, you can run 

3dMVM_validator -dataTable ILL.txt

With any luck, this will open up an AFNI window, and a page in your browser. Now, follow these instructions. You should be able to follow the instructions and make yourself a 3dMVM script. You should even be able to make one that works. ;-)

There are a few things I've noticed when playing with this. (Note, it's a pretty new feature, so these may be resolved by the time you read this.)

  • When I want to do a quantitative variable (for instance, does the TaskA_vs_TaskB contrast vary with Age?), 3dMVM_validator has me set this up with a "Age : 1" contrast code, but what I really want is "Age : ".

After running this, I find that I would often like to have the FDR correction in the file (but maybe that's just me.) To add that, run 

3drefit -addFDR MVM_output.nii

with "MVM_output.nii" being, of course, the output of your 3dMVM script. Now, you should be able to view the statistics with FDR correction.

ICA using melodic

To do an ICA, we're going to use the FSL tool melodic. We're going to give melodic the data that we cleaned up earlier using AFNI. You can use melodic on the command line or the gui. First we'll walk through the gui.

In a terminal window (assuming you've installed FSL), 

 fsl

This will bring up the FSL GUI, and you can choose MELODIC ICA. When this GUI comes up, we'll start with the "Select 4D data," and have it point to the cleaned up data we want it to use. This will be the errts dataset we generated earlier with just a few tweaks. First, AFNI outputs in .BRIK and .HEAD files, and we need a .nii file for FSL to read. Also, one of the first things melodic is going to want to do is determine what's brain and what's not. The easy way to do that is to look for high-intensity voxels near the center of the image. BUT, we already de-meaned our data, so all the voxels have a mean of zero. To "fix" that for melodic, we're simply going to add some height to the brain voxels. For this, you're going to need a brain mask, which should have been generated when we did the afni_proc.py steps earlier. We'll use mask_epi_anat.*+tlrc 

 (in the *results folder)
 3dcalc -a mask_epi_anat.${sid}+tlrc -b errts.${sid}.fanaticor+tlrc -expr 'b+100*a' -prefix tall.$sid.nii -datum float

In this step, we added 100 to all the in-brain areas, and also converted to .nii. (The -datum float bit is because the mask is type byte, but the errts is type float, and we need our output to be type float. AFNI will probably guess right, but let's make sure.) So now, we're ready to roll.

In the Melodic window, select "Select 4D data" and navigate to tall.$sid.nii. Once you've selected your data, it will automatically detect the TR and the number of volumes. Make sure that looks right. Now, we'll move on to the next tab, which is "Pre-Stats".

We've already processed this, so now we're just going to tell Melodic a big 'ol Thanks, but No Thanks. So, we're setting motion correction to "None," Slice timing correction to "None," BET brain extraction (masking) to "None," Smoothing to "0," and turning off Highpass filtering. (I know, I know. If we were just going to turn BET off, why'd we bother with the mask thing, right? Well, it's still going to find the brain bits, it's just not going to use BET to do it... or, more accurately, I believe it's still going to set all the nearly-zero voxels to be zero.)

On the next page, turn registration off.

Here, it's a bit trickier. Both of these options are reasonable options.

Variance normalise [sic] : I almost alwasy turn this off. The idea here is that some parts of the brain are brighter and others are dimmer just because of where they are in relation to the coils or the sinuses or whatever. In your raw data, variance is pretty much proportional to signal, so when you variance-normalize, you are giving all the voxels in the brain equal votes, regardless of how bright they are. Here, we no longer have a raw signal, rather, we have a residual timecourse from the nuisance regressors. I'm inclined to say that if there isn't so much data left that wasn't explained best by your noise regressors, let it be weighted less heavily in the ICA. If for your project, you're particularly concerned about areas that may have low signal, you could turn this on, but you could also calculate your residuals from the percent signal change in the afni_proc.py step. I think the latter is preferable.
Automatic dimensionality estimation : Here, you can choose how many components you want, or you can let Melodic choose for you based on the variance explained. Generally, I start with automatic detection, look at the output plot, and may run it again choosing the dimensions myself later. Normally, it's fine to do the automatic. But sometimes, you end up with 400 components or something, and it's just not so useful anymore.

The final tab is "Post-Stats." Here, we'll keep the Threshold IC maps checked, but set the threshold to 0.0. (We want it to give us the maps, but show us all of the output.) Set the Background image to "Mean functional" (because that's the only image you're giving it, even though it's not going to have a whole lot of spatial information). And check the Output full stats folder. Now, hit "Go."

The next thing you'll see is a webpage pop up with a Progress Report. This is really useful because it's printing out the commands it is using to run the analysis. These can now be borrowed for batch scripting later. Specifically, you care about the line that starts with /Applications/fsl/bin/melodic, because that's your melodic output.

More to come

Retrieved from "https://openwetware.org/mediawiki/index.php?title=Beauchamp:RestingState&oldid=1045216"