Regarding the overall pre-processing steps

AFNI Message Board
1

Dear AFNI experts,

Hello, I’m trying to determine the pre-processing steps for my data.
However, it is still not clear for me whether I’m right or not.
Also, there are several steps that I want to ask.

The below is the questions which I want to ask.

  1. [detrend part] I tried to conduct detrend after ‘scale’.
    However, when I check the ‘TSNR average’, the value seems too small.
    Thus, I also tried to conduct ‘detrend’ before the ‘scale’, and added the mean values from the ‘blur’ (the previous step result).
    That time, I could get the higher value of the ‘TSNR average’, but I their seems no big difference between those two method.
    Question: Is it the better way to increase ‘TSNR average’? If so, is it okay to conduct ‘detrend’ before ‘scale’?
    => related code (detrend before scale):

3dTstat -mask mask_group+tlrc -prefix mean.$subj.r$run.pb04 -mean pb04.$subj.r$run.blur+tlrc

# ================================ detrend =================================
3dDetrend -prefix pb05.$subj.r$run.detrend -polort 3 pb04.$subj.r$run.blur+tlrc

3dcalc -a pb05.$subj.r$run.detrend+tlrc -b mean.$subj.r$run.pb04+tlrc -c mask_group+tlrc -d full_mask.$subj+tlrc \
           -expr 'a*step(c)*step(d)+b*step(c)*step(d)' -prefix pb05.$subj.r$run.detrend_msk

# ================================= scale ==================================
# scale each voxel time series to have a mean of 100
# (be sure no negatives creep in)
# (subject to a range of [0,200])
foreach run ( $runs )
    3dTstat -prefix rm.mean_r$run pb05.$subj.r$run.detrend_msk+tlrc
    3dcalc -a pb05.$subj.r$run.detrend_msk+tlrc -b rm.mean_r$run+tlrc \
           -expr 'min(200, a/b*100)'       \
           -prefix pb06.$subj.r$run.scale
end

*2. [normalization related]
- I just curios how can I match the brain map to the MNI space more clearly.
The reason is, I use ‘@auto_tlrc -base MNI_avg152T1+tlrc -input anat_ns+orig -no_ss’.
After conducting @auto_tlrc and volreg part, I could get the normalized data.
However, when I check the mask (e.g. mask_group, MNI data), the individual map isn’t match perfectly.
I know that it is difficult to match perfectly, and little bit difference is not that severe issue.
but, in the ‘scale’ part, the default code usually use ‘mask_group’, and sometimes it seems larger than the individual brain.
Thus, it seems that the values which are non-brain part was included.
(I want to use the mask from the ‘mask_group’, instead of the individual mask, to get the similar brain from every individual subjects.
Also for the group analysis)
Those value could be meaningless one when we conduct the group analysis, but I want to know if there are some options,
that we could match the individual map more correctly into the template (e.g. MNI atlas).
Question: is there some options/functions that I could use for matching the individual brain into MNI atlas or mask_group?

The bottom will show the overall code which I use for the pre-processing (start from the despiking).
In advance, thank you for reading my questions and code.
It would be really appreciate if you give me some comments/advices.
Thank you.


# ================================ despike =================================
# apply 3dDespike to each run
foreach run ( $runs )
    3dDespike -NEW -nomask -prefix pb01.$subj.r$run.despike \
        pb00.$subj.r$run.tcat+orig
end
# ================================ retroicor =================================
# apply 3dDespike to each run
foreach run ( $runs )
    3dretroicor -prefix pb02.$subj.r$run.retroicor -resp /data/path/resp.1D -card /data/path/card.1D -order 2 \
        pb01.$subj.r$run.despike+orig
end
# ================================= tshift =================================
# time shift data so all slice timing is the same 
foreach run ( $runs )
    3dTshift -TR 2s -tzero 0 -slice 35 -quintic -tpattern alt+z2 -prefix pb03.$subj.r$run.tshift \
             pb02.$subj.r$run.retroicor+orig
end

# --------------------------------
# extract volreg registration base
3dbucket -prefix vr_base_min_outlier                           \
    pb03.$subj.r$minoutrun.tshift+orig"[$minouttr]"

# ================================= align ==================================
# for e2a: compute anat alignment transformation to EPI registration base
# (new anat will be intermediate, stripped, anat_ns+orig)
align_epi_anat.py -anat2epi -anat anat+orig      \
       -save_skullstrip -suffix _al_junk         \
       -epi vr_base_min_outlier+orig -epi_base 0 \
       -epi_strip 3dAutomask                     \
       -volreg_method 3dWarpDrive -giant_move	 \
       -volreg off -tshift off

# ================================== tlrc ==================================
# warp anatomy to standard space
@auto_tlrc -base MNI_avg152T1+tlrc -input anat_ns+orig -no_ss -dxyz 3.0 -OK_maxite -inweight -init_xform AUTO_CENTER
# I used -inweight -init_xform AUTO_CENTER for only one subject, because the brain was scaled and tilted strangely

# store forward transformation matrix in a text file
cat_matvec anat_ns+tlrc::WARP_DATA -I > warp.anat.Xat.1D

# ================================= volreg =================================
# align each dset to base volume, align to anat, warp to tlrc space

# verify that we have a +tlrc warp dataset
if ( ! -f anat_ns+tlrc.HEAD ) then
    echo "** missing +tlrc warp dataset: anat_ns+tlrc.HEAD" 
    exit
endif

# register and warp
foreach run ( $runs )
    # register each volume to the base
    3dvolreg -verbose -zpad 1 -base vr_base_min_outlier+orig    \
             -1Dfile dfile.r$run.1D -prefix rm.epi.volreg.r$run \
             -cubic                                             \
             -1Dmatrix_save mat.r$run.vr.aff12.1D               \
             pb03.$subj.r$run.tshift+orig

    # create an all-1 dataset to mask the extents of the warp
    3dcalc -overwrite -a pb03.$subj.r$run.tshift+orig -expr 1   \
           -prefix rm.epi.all1

    # catenate volreg/epi2anat/tlrc xforms
    cat_matvec -ONELINE                                         \
               anat_ns+tlrc::WARP_DATA -I                       \
               anat_al_junk_mat.aff12.1D -I                     \
               mat.r$run.vr.aff12.1D > mat.r$run.warp.aff12.1D

    # apply catenated xform: volreg/epi2anat/tlrc
    3dAllineate -base anat_ns+tlrc                              \
                -input pb03.$subj.r$run.tshift+orig             \
                -1Dmatrix_apply mat.r$run.warp.aff12.1D         \
                -mast_dxyz 3                                    \
                -prefix rm.epi.nomask.r$run			\
		-lpc -automask -source_automask

    # warp the all-1 dataset for extents masking 
    3dAllineate -base anat_ns+tlrc                              \
                -input rm.epi.all1+orig                         \
                -1Dmatrix_apply mat.r$run.warp.aff12.1D         \
                -mast_dxyz 3 -final NN -quiet                   \
                -prefix rm.epi.1.r$run			\
		-lpc -automask -source_automask

    # make an extents intersection mask of this run
    3dTstat -min -prefix rm.epi.min.r$run rm.epi.1.r$run+tlrc
end

# make a single file of registration params
cat dfile.r*.1D > dfile_rall.1D

# ----------------------------------------
# create the extents mask: mask_epi_extents+tlrc
# (this is a mask of voxels that have valid data at every TR)
# (only 1 run, so just use 3dcopy to keep naming straight)
3dcopy rm.epi.min.r01+tlrc mask_epi_extents

# and apply the extents mask to the EPI data 
# (delete any time series with missing data)
foreach run ( $runs )
    3dcalc -a rm.epi.nomask.r$run+tlrc -b mask_epi_extents+tlrc \
           -expr 'a*b' -prefix pb04.$subj.r$run.volreg
end

# warp the volreg base EPI dataset to make a final version
cat_matvec -ONELINE                                             \
           anat_ns+tlrc::WARP_DATA -I                           \
           anat_al_junk_mat.aff12.1D -I  > mat.basewarp.aff12.1D

3dAllineate -base anat_ns+tlrc                                  \
            -input vr_base_min_outlier+orig                     \
            -1Dmatrix_apply mat.basewarp.aff12.1D               \
            -mast_dxyz 3                                        \
            -prefix final_epi_vr_base_min_outlier		\
	    -lpc -automask -source_automask

# create an anat_final dataset, aligned with stats
3dcopy anat_ns+tlrc anat_final.$subj

# record final registration costs
3dAllineate -base final_epi_vr_base_min_outlier+tlrc -allcostX  \
            -input anat_final.$subj+tlrc -lpc -automask -source_automask |& tee out.allcostX.txt

# -----------------------------------------
# warp anat follower datasets (affine)
3dAllineate -source anat+orig                                   \
            -master anat_final.$subj+tlrc                       \
            -final wsinc5 -1Dmatrix_apply warp.anat.Xat.1D      \
            -prefix anat_w_skull_warped				\
	    -lpc -automask -source_automask

# ================================== blur ==================================
# blur each volume of each run
foreach run ( $runs )
    3dBlurToFWHM -FWHM 8 -automask				\
            -ACF -rate 0.2 -temper				\
            -input pb04.$subj.r$run.volreg+tlrc			\
            -prefix pb05.$subj.r$run.blur
end

# ================================== mask ==================================
# create 'full_mask' dataset (union mask)
foreach run ( $runs )
    3dAutomask -dilate 1 -prefix rm.mask_r$run pb05.$subj.r$run.blur+tlrc
end

# create union of inputs, output type is byte
3dmask_tool -inputs rm.mask_r*+tlrc.HEAD -union -prefix full_mask.$subj

# ---- create subject anatomy mask, mask_anat.$subj+tlrc ----
#      (resampled from tlrc anat)
3dresample -master full_mask.$subj+tlrc -input anat_ns+tlrc           \
           -prefix rm.resam.anat

# convert to binary anat mask; fill gaps and holes
3dmask_tool -frac 1.0 -fill_holes -input rm.resam.anat+tlrc  \
            -prefix mask_anat.$subj

# compute overlaps between anat and EPI masks
3dABoverlap -no_automask full_mask.$subj+tlrc mask_anat.$subj+tlrc    \
            |& tee out.mask_ae_overlap.txt

# note Dice coefficient of masks, as well
3ddot -dodice full_mask.$subj+tlrc mask_anat.$subj+tlrc               \
      |& tee out.mask_ae_dice.txt

# ---- create group anatomy mask, mask_group+tlrc ----
#      (resampled from tlrc base anat, MNI_avg152T1+tlrc)
3dresample -master full_mask.$subj+tlrc -prefix ./rm.resam.group      \
           -input /home/hdw/abin/MNI_avg152T1+tlrc

# convert to binary group mask; fill gaps and holes
3dmask_tool -frac 1.0 -fill_holes -input rm.resam.group+tlrc \
            -prefix mask_group

3dresample -master anat_final.$subj+tlrc -prefix new_mni -inset MNI_avg152T1+tlrc
3dSeg -anat anat_final.$subj+tlrc -mask new_mni+tlrc

3dcopy Segsy/Classes+tlrc .

# ================================= scale ==================================
# scale each voxel time series to have a mean of 100
# (be sure no negatives creep in)
# (subject to a range of [0,200])
foreach run ( $runs )
    3dTstat -prefix rm.mean_r$run pb05.$subj.r$run.blur+tlrc
    3dcalc -a pb05.$subj.r$run.blur+tlrc -b rm.mean_r$run+tlrc \
           -c mask_group+tlrc                                  \
           -expr 'c * min(200, a/b*100)*step(a)*step(b)'       \
           -prefix pb06.$subj.r$run.scale
end

# ================================ detrend =================================
3dDetrend -prefix pb07.$subj.r$run.detrend -polort 3 	\
	pb06.$subj.r$run.scale+tlrc

3dcalc -a Classes+tlrc -expr 'equals(a,1)' -prefix CSF
3dcalc -b Classes+tlrc -expr 'equals(b,3)' -prefix WM

3dresample -master pb07.$subj.r$run.detrend+tlrc -inset WM+tlrc -prefix WM_Resampled
3dresample -master pb07.$subj.r$run.detrend+tlrc -inset CSF+tlrc -prefix CSF_Resampled
3dmaskave -quiet -mask WM_Resampled+tlrc pb07.$subj.r$run.detrend+tlrc > WM_Timecourse.1D
3dmaskave -quiet -mask CSF_Resampled+tlrc pb07.$subj.r$run.detrend+tlrc > CSF_Timecourse.1D
2

Hello,

  1. If the mean is left in, detrending should be irrelevant to
    the TSNR, or even to the regression.

TSNR is computed after the regression (in afni_proc.py), which
applies trend parameters. Detrending should also not affect the
regression betas or the statistics, as far as I can see, since
all statistics are computed above and beyond the baseline.

One way to improve TSNR is to get your subjects to move less. :slight_smile:

  1. The group_mask dataset is used for nothing, by default. It
    is not a great mask to use in general.

The mask that is applied in scaling is the extents mask, which
should only affect a limited number of voxels.

So yes, you should see data outside the brain. We are not
intending to mask out data at the single subject level. That
can be saved for whatever you use as a group mask.

Do not evaluate registration by looking at masks, but by looking
at the actual datasets in question (anat/EPI at single subject
level, and anat/template at group level).

  • rick
3

Dear Rick,

Thank you for the detailed comments/advices.
It really helped to understand more about the TSNR and the mask related part.
Again, thank you so much.

Regards,
DaWoon

4

Dear AFNI expert,

I have three additional question related to the tshift, @auto_tlrc, and wm/csf regressing part.
(My apologies in advance, it could be easy/irritating questions…:-()
The reason is, I thought that it is okay to use, but I heard it is difficult to say correct or not to use in pre-processing from my supervisor…
Thus, may I ask you whether it is okay to use?
The code is shown below…

  • [tshift] I defined the TR and the number of slice and the sequence that I acquired during scanning.
    . (Actually, I heard that we don’t need to set TR, the number of slices, etc… but I just wanted to define clearly…)

# ================================= tshift =================================
# time shift data so all slice timing is the same 
foreach run ( $runs )
    3dTshift -TR 2s -tzero 0 -slice 35 -quintic -tpattern alt+z2 -prefix pb03.$subj.r$run.tshift \
             pb02.$subj.r$run.retroicor+orig
end

-[tlrc] there are two ways, and the same thing is that I just set the dimension of xyz…

  1. for the odd subject which brain scaled differently after tlrc, and the other one is for the normal subject

# ================================== tlrc ==================================
# warp anatomy to standard space
@auto_tlrc -base MNI_avg152T1+tlrc -input anat_ns+orig -no_ss -dxyz 3.0 -OK_maxite -inweight -init_xform AUTO_CENTER
  1. for the normal subject

# ================================== tlrc ==================================
# warp anatomy to standard space
@auto_tlrc -base MNI_avg152T1+tlrc -input anat_ns+orig -no_ss -dxyz 3.0 -OK_maxite
  • Extracting signals from the WM/CSF to use in the regressors

# To match the dimension..
3dresample -master anat_final.$subj+tlrc -prefix new_mni -inset MNI_avg152T1+tlrc
3dSeg -anat anat_final.$subj+tlrc -mask new_mni+tlrc

3dcopy Segsy/Classes+tlrc .

3dcalc -a Classes+tlrc -expr 'equals(a,1)' -prefix CSF
3dcalc -b Classes+tlrc -expr 'equals(b,3)' -prefix WM

3dresample -master pb07.$subj.r$run.detrend+tlrc -inset WM+tlrc -prefix WM_Resampled
3dresample -master pb07.$subj.r$run.detrend+tlrc -inset CSF+tlrc -prefix CSF_Resampled
3dmaskave -quiet -mask WM_Resampled+tlrc pb07.$subj.r$run.detrend+tlrc > WM_Timecourse.1D
3dmaskave -quiet -mask CSF_Resampled+tlrc pb07.$subj.r$run.detrend+tlrc > CSF_Timecourse.1D

. After conducting this one, I input WM_Timecourse.1D and CSF_Timecourse.1D into the GLM regress part.

Thank you for reading and thank you in advance…

Regards,
DaWoon

5

Hi DaWoon,

The tshift step is fine. Many people do that, many do not.
For shorter events, one might hope to get better temporal
alignment with the BOLD signal that way.

It is strange to have both -tzero 0 and -slice 35 in there.
Only one option specifying what time to align to is applied.

Addintion options for @auto_tlrc to succeed is fine. The
AUTO_CENTER option is generally only needed if the coordinates
in the dataset are not accurate. Highlighting -dxyz 3.0
suggests that the voxel size of the odd subject is different.
If so, and if you can afford dropping a subject, it might be
good to drop that one.

The CSF mask from 3dSeg is not good to use on its own, as it
includes CSF from all over the brain (not just ventricles).
Consider Example 11b from “afni_proc.py -help”. It passes
in a ventricle mask from the group template, and intersects
it with the 3dSeg CSF mask to make a subject ventricle mask.

  • rick