静息态的数据处理现在都比较简单了,因为有了fmriprep和qsiprep流程工具,没有特别的难度。
而关于任务态,则有一些独特的处理方式,因为最近要做任务态的数据分析,所以学习一下。
先学习下R的4个重输出函数:cat , sink, wirteLines, write.table
1.cat函数即能输出到屏幕,也能输出到文件.
使用方式:cat(… , file = “”, sep = " ", fill = FALSE, labels = NULL,append = FALSE)
有file时,输出到file。无file时,输出到屏幕。
append参数:布尔值。TRUE时,输出内容追加到文件尾部。FALSE,覆盖文件原始内容。
2.sink函数将输出结果重定向到文件。
使用方式:sink(file = NULL, append = FALSE, type = c(“output”, “message”),split = FALSE)
append参数:布尔值。TRUE时,输出内容追加到文件尾部。FALSE,覆盖文件原始内容。
3.writeLines函数将字符串向量输出到文件中(会覆盖原始内容)
使用方式:writeLines(text, con = stdout(), sep = “\n”, useBytes = FALSE)
text:字符串向量;con:输出文件
4.write.table()函数将dataframe的内容输出到文件。
事件
静息态一般不考虑事件,而任务态有多个block的事件设计,先把这些block合并成1个bold文件,如果不会用代码可以用fslmath -t的命令。
预处理
等bold图像按照常规fmriprep处理以后,这个就是预处理完的准备好的数据,用于后续统计分析和指标计算。
统计建模
任务态的一个特点就是事件相关,需要有events表格文件,注明每个时间点属于哪个事件。
拟合函数:
The Hemodynamic Response Function (HRF)
在患者被刺激后,其bold信号一般经历“逐步增加-至顶峰-再降低”的过程,这个过程的形状可以用gamma分布的函数去拟合,
伽玛分布(Gamma Distribution)是统计学的一种连续概率函数,是概率统计中一种非常重要的分布。“指数分布”和“χ2分布”都是伽马分布的特例。
当使用 Gamma 分布去拟合 BOLD 响应时,可将其称为典型血液动力学响应函数或HRF,而gamma分布又可以称为基函数(basis function)。
stick function
gamma分布去拟合的bold信号,其刺激一般是数秒内的。如果是持续了一段时间的刺激,如持续刺激15秒,或者持续完成某个事件15秒以上,这个时候可以用另一个函数去拟合,stick function,它是gamma和block刺激进行卷积的结果。
上图中的红色和绿色,取其平均值,这就是一级分析,SPM里面的first-level analysis。
一级分析就是对每个被试进行建模,二级分析是利用一级分析产生的建模参数,在GLM上进行平均和对比估计,contrast。
不过不同的是,SPM的second-level是组分析,而fsl的higher-level才是组分析,fsl的second-level只是对估计的参数进行平均。
GLM模型参数选择:
Fixed Effects:
Mixed Effects: Simple OLS (Ordinary Least Squares): This will perform a t-test on the average parameter estimates calculated for each subject, without taking into account the variability between the runs for each subject;
Mixed Effects: FLAME 1: Weight each subject’s parameter estimate by the variance of that contrast estimate. In other words, a subject with relatively low variance will be weighted more, and a subject with relatively high variance will be weighted less;
Mixed Effects: FLAME 1+2: A more rigorous version of FLAME 1. It takes much longer, and is only helpful for analyzing small samples (e.g., 10 subjects or fewer);
Randomise: A non-parametric test.
使用混合效应模型,去拟合GLM模型,然后对voxel或者cluster进行校正,就得到最后的结果。即事件刺激的脑区。
Finite Impulse Response (FIR)
有限脉冲响应
HRF是最常用的拟合事件的函数,它的特点是为特定事件的所有被试估计单个最佳拟合参数。
假设需要观察事件开始后bold信号随活动的时间进程而变化,就需要使用FIR函数去拟合。
FIR函数,有限脉冲响应模型,我们可以指定时间窗口的长度和要估计的时间点数,例如在一个事件从开始到持续到完成以后的10秒内,我们想观察这10秒的每2秒不同阶段变化,可以使用FIR函数生成5个拟合参数,而HRF对于这10秒只能生成1个拟合参数。其实可以类比于动态分析里的滑动时间窗。
(Prefrontal cortex organization: dissociating effects of temporal abstraction, relational abstraction, and integration with FMRI, Cerebral cortex, 2013.)
如上图,Y轴是beta估计值,代表激活,X轴是分隔的时间段(2秒)。
(Exploring fMRI Results Space: 31 Variants of an fMRI Analysis in AFNI, FSL, and SPM, 2016)
Nilearn完成任务态分析
Nilearn的相关函数:
Nilearn的一级分析和二级分析,跟SPM是对应的,一级在被试上,二级是组分析。
# fit GLM
first_level_model = FirstLevelModel(t_r=t_r, slice_time_ref=slice_time_ref)
first_level_model = first_level_model.fit(
run_imgs=adhd_dataset.func[0], design_matrices=design_matrix
)
# contrast estimate
z_map = first_level_model.compute_contrast(
contrasts["seed_based_glm"], output_type="z_score"
)
# generate report
from nilearn.reporting import make_glm_report
report = make_glm_report(
first_level_model,
contrasts=contrasts,
title="ADHD DMN Report",
cluster_threshold=15,
min_distance=8.0,
plot_type="glass",
)
nilearn.glm.first_level.FirstLevelModel(t_r=None, slice_time_ref=0.0, hrf_model=‘glover’, drift_model=‘cosine’, high_pass=0.01, drift_order=1, fir_delays=[0], min_onset=-24, mask_img=None, target_affine=None, target_shape=None, smoothing_fwhm=None, memory=Memory(location=None), memory_level=1, standardize=False, signal_scaling=0, noise_model=‘ar1’, verbose=0, n_jobs=1, minimize_memory=True, subject_label=None, random_state=None)
hrf_model: ‘spm’ (default), ‘spm + derivative’, ‘spm + derivative + dispersion’, ‘glover’…‘fir’
在大多数情况下,spm标准模型和glover模型不会产生大的差异。
固定效应统计:compute_fixed_effects
设计矩阵:
nilearn的二级分析:
用make_second_level_design_matrix函数。
经过一级二级分析,然后再进行校正,就可以得到全脑激活图了。
使用nilearn的另一个好处就是,除了统计分析,还可以直接用这些结果进行机器学习。“用于基于任务的功能连接和解码的 Beta 系列建模”