fMRI & R

R语言的统计功能如此强大，我在想能不能用R来分析fMRI数据，于是google了一下。

有一个关于用R分析fMRI数据的课程：Neuroimaging Analysis within R

R有一个叫做“fmri”的package。（具体怎么用，以后再说）

R还有一个叫做“neurobase”的package。

library(neurobase)
setwd("~/Desktop/validating-fmri/data")
data <- readnii("sub-10159/sub-10159_task-rest_bold_space—MNI152NLin2009cAsym_preproc.nii.gz")
#The function readnii requires one argument: the name of the file we want to read. 
#What comes out is an object of type oro.nifti. 
#Essentially, it is a 4D array, with extra metadata specific to fMRI data. 

#We can therefore apply functions that we can apply to arrays. Let’s look for example at the shape of our data.
dim(data)
[1]  65  77  49 152
#This shows us that the first three dimensions are 65 x 77 x 49. 
#These correspond to the x-, y- and z-coordinates of the brain. 
#The fourth dimension is time, which means there are 152 time points measured.

#We can also run a function from neurobase that only works on this type of object:
check_nifti(data)
NIfTI-1 format
  Type            : nifti
  Data Type       : 16 (FLOAT32)
  Bits per Pixel  : 32
  Slice Code      : 0 (Unknown)
  Intent Code     : 0 (None)
  Qform Code      : 2 (Aligned_Anat)
  Sform Code      : 1 (Scanner_Anat)
  Dimension       : 65 x 77 x 49 x 152
  Pixel Dimension : 3 x 3 x 4 x 2
  Voxel Units     : Unknown
  Time Units      : sec
#This shows some information that is stored in the “header” of the file. For example, it tells us the pixel dimensions. 
#The first three numbers tell us the size of the voxels in space: 3 x 3 x 4 cm.
#The last number tells us the dimension in time: a scan was taken every 2 seconds.

#Another interesting function is the visualisation of nifti’s. 
#Note that we only visualise one timepoint, the first one. R will automatically show the first timepoint.
orthographic(data)
orthographic(data,xyz=c(20,20,30))
#The argument xyz allows us to specify where to draw the cross-hairs

#To look at the value of one specific point in time and space, we can use indexes.
data[20,20,30,1]
[1] 1148.524
#（20,20,30表示坐标，1表示时间点）

#What is more interesting is to look at how the value changes over time, which we can do by omitting the index in the 4th dimension:

data[20,20,30,]
[1] 1148.524 1148.461 1142.563 1141.060 1127.099 1131.675 1135.632 1128.941
  [9] 1160.061 1161.512 1137.481 1116.822 1109.309 1109.323 1093.169 1110.955
 [17] 1109.561 1116.910 1134.975 1118.357 1132.719 1150.193 1165.766 1173.594
 [25] 1176.839 1142.286 1145.902 1159.892 1152.587 1136.450 1127.974 1148.936
 [33] 1152.663 1174.247 1168.272 1147.255 1159.906 1184.832 1145.608 1141.148
 [41] 1132.086 1137.863 1133.563 1143.711 1141.173 1128.889 1124.897 1132.619
 [49] 1144.724 1153.767 1133.262 1135.230 1137.793 1165.782 1160.442 1159.418
 [57] 1165.768 1156.091 1169.558 1182.266 1167.106 1157.483 1163.735 1131.992
 [65] 1121.514 1111.989 1121.143 1141.684 1122.967 1126.063 1121.025 1138.645
 [73] 1132.201 1133.171 1132.688 1134.140 1127.046 1127.342 1134.374 1148.151
 [81] 1147.634 1145.883 1143.286 1142.210 1138.209 1123.161 1148.567 1149.430
 [89] 1155.253 1115.818 1116.521 1119.263 1136.349 1144.674 1151.027 1158.472
 [97] 1145.075 1124.284 1120.904 1129.625 1131.517 1120.666 1126.623 1109.636
[105] 1093.798 1126.940 1126.045 1120.136 1133.749 1101.970 1098.026 1116.716
[113] 1142.895 1143.010 1138.357 1127.057 1102.292 1122.499 1136.995 1126.363
[121] 1137.567 1122.647 1142.551 1131.645 1125.780 1120.023 1122.369 1132.125
[129] 1104.158 1107.103 1088.539 1116.694 1149.974 1146.161 1140.482 1114.848
[137] 1099.089 1121.268 1141.496 1141.877 1129.998 1106.165 1107.426 1105.125
[145] 1089.509 1103.882 1114.204 1113.552 1096.908 1104.983 1102.881 1121.604
#If we plot this vector, we can see how the measured signal in this voxel changes over time.
plot(data[20,20,30,],type='l')