Matilde QC

Matilde summarized the main MEICA quality checks implement following Kundu et al.:

Please download this Word document for a comprehensive documentation.

Quality check after meica preprocessing

Meica’s output 1

medn.nii.gz 

'Denoised' BOLD time series after basic preprocessing, T2* weighted averaging and ICA denoising.

Meica’s output 2

tsoc.nii.gz 

'Raw' BOLD time series dataset after basic preprocessing and T2* weighted averaging of echoes (i.e. 'optimal combination'). 'Standard' denoising can be assessed on this dataset (e.g. motion regression, physio correction, scrubbing) for comparison to ME-ICA denoising.

Meica’s output 3

mefc.nii.gz 

Component maps (in units of $\delta S$) of accepted BOLD ICA components. Use this dataset for ME-ICR seed-based connectivity analysis.

Coregistration

This step verifies that the mprage produced after preprocessing is aligned to the functional timeseries. Select the following images to check successful coregistration.

mprage_ns_atnl.nii.gz

output file produced by meica corresponding to skullstripped anatomical (ns) to which non linear warp (nl) to standard space (at) is applied

resting_medn.nii.gz

denoised time series after basic preprocessing

By automatizing this script, it will be possible to generate visualization of the overlap for each subject.

afni -com 'OPEN_WINDOW A.sagittalimage' -com "SWITCH_UNDERLAY mprage_ns_atnl.nii.gz" -com "SWITCH_OVERLAY ${sub}resting_medn.nii.gz" -com "SAVE_JPEG A.sagittalimage Coreg_${sub}.jpg" -com "QUIT"

It is necessary to verify that the two images overlap. If this is not the case for some subjects, solutions need to be identified. A first easy try entails rerunning the individual subjects for which overlap is not accomplished. For example, the image below shows a good co-registration for 10 subjects.

Inspect number and quality of BOLD-like components identified

This step allows visual inspection of the components identified by meica and included in the denoised medn file. In general, as a rule of thumb, a subject is good if >10 components are identified (Check with Petra, Kundu P. gave this number I believe at the very early stages). Lower number indicates that many artefacts might have contaminated the data of the subject to begin with. Exclusion of the participants with less than 10 components should be considered.

Number of components retained for each subject can be simply verified by inspecting the ctab.txt file generated by meica and included in subject’s folder. The number of accepted BOLD components is indicated as:

#No. accepted BOLD-like components, i.e. effective degrees of freedom for correlation (lower bound; DFe)

With this step, the user will also be able to visualize the components identified and retained as BOLD signal. Those should roughly correspond to known functional networks. Example of the components are shown in Kundu et al., PNAS 2013, Figure S3. The components are ordered in descending order from largest to smallest amount of variance explained.

Some steps are required to visualize those components

• Apply standard space transformation to relevant file using the nifti_tool command which can be found at the end of the .sh file generated by meica (Kundu P, personal communication).

The relevant file is feats_OC2.nii (I think this file indexes the components not transformed in standard space but some better explanation should be provided) and it can be found in the TED subfolder of each subject. (Safer to first copy the relevant file to a separate directory, apply the transformation on the copied file, and check components on that file). The following transformation should be applied (previous verification that this nifti_tool is the same that can be found at the end of .sh file generated in the directory of each subject)

nifti_tool -mod_hdr -mod_field sform_code 2 -mod_field qform_code 2 -infiles ${sub}_feats_OC2.nii –overwrite

• Create txt file with the number of accepted components for each subject (this command counts the number of commas, so factually this file will include number of components – 1, see step c). File accepted.txt is found in TED subfolder of each subject

fgrep -o , accepted.txt | wc -l > ${sub}_ncomp.txt

• Create output for each component. By using 3dclust command, a table is created reporting all x y z peak coordinates of significant activation at a given threshold, for all the accepted components. (Seq 0 $n, is implemented to account for all components from step b)

n=`cat *ncomp.txt`
for i in $(seq 0 $n)
do
3dclust -1Dformat -nosum -1dindex ${i} -1tindex ${i} -2thresh -3.0 3.0 -dxyz=1 1.01 20 ${sub}_feats_OC2.nii >ReportTableComp_${sub}_${i}.txt
done

• Extract coordinates for main peak activation for each component (i) of every subject (sub)

n=`cat *ncomp.txt`
for i in $(seq 0 $n)
do 
sed -n "13p" ReportTableComp_${sub}_${i}.txt >RowCoord_${sub}_${i}.txt
done

• Generate figure for each component at peak coordinates x y z and display it with AFNI command at a given standard threshold

n=`cat *ncomp.txt`
for i in $(seq 0 $n)
do 
awk '{print $14 }' RowCoord_${sub}_${i}.txt > x.txt			
awk '{print $15 }' RowCoord_${sub}_${i}.txt > y.txt
awk '{print $16 }' RowCoord_${sub}_${i}.txt > z.txt

x=`cat *x.txt`
y=`cat *y.txt`
z=`cat *z.txt`

afni -com 'OPEN_WINDOW A.axialimage' -com "SWITCH_UNDERLAY ${sub}_mprage_ns_atnl.nii.gz" -com "SWITCH_OVERLAY ${sub}_feats_OC2.nii" -com "SET_SUBBRICKS A 0 ${i} ${i}" -com 'SET_THRESHNEW A 3.0 *p' -com "SET_DICOM_XYZ A $x $y $z " -com "SAVE_JPEG A.axialimage SnapComp_${sub}_${i}.jpg" -com "QUIT"
done

Display rigid-body movement parameters traces

See previous section on FD for recommendation on dealing with excessive motion. Example displaying rigid-body movement parameters can be found in Kundu et al., PNAS 2013 Figure 1 top row.

1dplot -volreg -plabel ${sub} -jpeg ${sub}_1Dplot ${sub}_motion.1D

Display DVARS for high-k and low-k components

Delta variation signal (DVARS), computed as the root mean square average of the first derivatives of fMRI percent signal change time courses, identifies time points with rapidly changing fMRI signal (Menon et al., Magn Reson Med 1993). Correspondence between FD and DVARS indicates contamination of fMRI time series with motion artefact. At this step, the user can compute DVARS for ME-ICA high K time series data (putatively indicating BOLD components) and for low k time series data (putatively indicating non-BOLD components) without explicit spectral filtering. Reduction in variance of the DVARS trace after denoising compared to before, as well as reduction in correspondence between DVARS and FD traces, is indicating of good denoising performance. See Kundu et al., PNAS 2013 Figure 1. These steps were implemented after personal correspondence with Kundu P, limited explanation was given but it made intuitively sense to me.

• Compute mean of each voxel time series for the tsoc file

3dTstat -mean -prefix ${sub}_echo1234_tsoc_mean.nii.gz ${sub}_echo1234_tsoc.nii.gz

• Add high k components (file found in TED folder) to mean tsoc (high k file does not have the mean image in it; it is necessary to compute the mean image of tsoc and add it to high k).

3dcalc -a hik_ts_OC.nii -b ${sub}_echo1234_tsoc_mean.nii.gz -expr 'a+b' -prefix hik_ts_OC_m.nii.gz

• The same should be done for the low k components

3dcalc -a lowk_ts_OC.nii -b ${sub}_echo1234_tsoc_mean.nii.gz -expr 'a+b' -prefix lowk_ts_OC_m.nii.gz

• Compute dvars on the newly calculated high k file including the mean of tsoc, low k file including the mean of tsoc, and for the original tsoc file

python /home/mmsv2/meica/dvars.py hik_ts_OC_m.nii.gz
python /home/mmsv2/meica/dvars.py lok_ts_OC_m.nii.gz
python /home/mmsv2/meica/dvars.py ${sub}_echo1234_tsoc.nii.gz

(Matilde: I think the dvars.py file was sent by Prantik and not sure it is currently included in meica. Also notice that from correspondence with Petra there might be differences in the way DVARS is calculated depending on whether scripts come from Ameera vs. Prantik)

Plot the obtained DVARS and plot to verify DVARS is lower for hik components similarly to Kundu et al., PNAS 2013, Figure 1 and uncoupled from FD traces in the case of hik components.