Two levels' ICA can disclose steady-state and transient profiles

The code below is dependent on GIFT, FSL and AFNI. It shows the 2nd level ICA since the 1st level can be done through GIFT, FSL melodic etc..

Transient information can be retrieved by using the saved subject's profile as weight (per the analysis of the paper) to be multiplied to the dynamic profiles given by GIG-ICA or any other back-reconstruction method (such as dual regression).

1. #!/usr/bin/env bash
2. #
3. # gicaS.sh
4. #% Init: 2013-11-06.13:37
5. #% Version: 20170614.17:28 Xiaowei Song
6. #% Copyright (C) 2013~2020 Xiaowei.Song <dawnwei.song@gmail.com>
7. #% http://restfmri.net/forum/XiaoweiSong
8. #% Distributed under terms of the AFL (Academy Free license).
9. #
10. CALLDIR="$( cd "$( dirname "${0}" )" && pwd )"
11. elog(){ echo "$@" 1>&2 ; }
12.  echo "$(date) | $0 $*" >> $(basename $0).LOG
13. source $(dawnbin)/slib.bash.sh
14. 
     
15. shopt -s expand_aliases; source $(dawnbin)/slib.stack.sh
16. usage(){ printf "
17. Usage: (%s)
18.     ${0} [options]
19.     options:
20.     Spacial ICA by using GIFT
21.     1. if multi nifti files input, default do Zconcat to make them in spatial dim
22.     2. need manually concat in T dim if all inputs are 3D nifti ( such as thickness ) files
23. " "$(grep -m1 '#% Version:' `readlink -f ${0}`)" 1>&2 ; }
24. while [ $# -gt 0 ]; do
25.     case "$1" in
26.         '-v') VERBOSE=1 ; shift ;;
27.         '-d') DEBUG=1 ; shift ;;
28.         '-m'|'--mask') MASK=$2 ; shift 2 ;;
29.         '-dm'|'--rm-mean') DEMEAN=$2; shift 2 ;;
30.         -k|'-n') COMPNUM=$2; shift 2 ;;
31.         #usually be 99->30
32.         -evp|--expected-var-percent) EXPECTED_VAR_PERCENT=$2; shift 2;;
33.         '-p'|'--prefix') PREFIX=$2; shift 2 ;;
34.          --mst) MST=$2; shift 2;;
35.         -rd|--cut-dir) CUT_DIR=$2; shift 2;;
36.         '-*') elog 'Unknown parameters'; usage; exit 0 ;;
37.         *) break ;;
38.     esac
39. done
40. verbose=${VERBOSE:-0};
41. 
     
42. if [ $# -eq 0 ] || ! IsFileExist $MASK ; then usage; exit 0; fi
43. nii=$1
44. bnii=$(basename $nii .nii)
45. #%parse ICA options
46. OPT4ICA=${ICAOPTIONS:-{}}
47. fpca="${PREFIX:-${bnii}}_pca.nii"
48. rm -f $fpca #clean, must do  since SPM8 will only overwrite certain time points
49. 
     
50. if ! IsFileLoadable $MASK; then
51.     elog "Mask not loadable: $MASK"
52.     #automask based on first input nii
53.     MASK=$(getmp).nii
54.     3dAutomask -prefix $MASK $1
55. fi
56. read ax ay az <<<$(3dinfo -ni -nj -nk $MASK)
57. 
     
58. 
     
59. mkStack inputs
60. mkStack ftmps
61. #build all input and zstack
62. lst=${PREFIX:-jointICA}.lst
63. while [[ $# -gt 0 ]]; do
64.     nii=$1
65.     pushs inputs $(rmExtDir4nii $nii)
66.     gz=${nii##*.}
67.     prefix=$(rmExtDir4nii $nii)
68.     if [[ gz == "$gz" ]]; then
69.         nii=$(getmp).nii
70.         pushs ftmps $nii
71.         gunzip -c $1 > $nii
72.     fi
73.     echo $(readlink -f $nii) ; shift;
74. done > $lst
75. showStack inputs
76. #exit
77. nnii=$(cat $lst|wc -l)
78. 
     
79. #build jICA input
80. nii=$(getmp).nii  ; pushs ftmps $nii
81. msk=$(getmp).nii ; pushs ftmps $msk
82. 3dhdr $(cat $lst |xargs) $MASK
83. if [[ 1 -eq $nnii ]]; then
84.     3dcopy $(cat $lst |xargs) $nii #thus support nii.gz if only 1 input
85.     3dcopy $MASK $msk
86. else
87.     #3dZcat needs at least 2 inputs
88.     3dZcat -prefix $nii $(cat $lst |xargs)
89.     3dZcat -prefix $msk $(yes $MASK |head -n $nnii |xargs)
90. fi
91. 
     
92. 
     
93. result=${PREFIX:-${bnii}}_component_ica_.nii
94. resultTC=${PREFIX:-${bnii}}_timecourses_ica_.txt
95. rm $result $resultTC
96. 
     
97. tmp=$(getmp).m
98. cat > $tmp <<eom
99. try
100. 
      
101. [msko, hinfo]=icatb_read_data('$msk');
102. msk=find(msko);
103. 
      
104. data=icatb_read_data('$nii',[],msk);
105. data=icatb_preproc_data(data, '${DEMEAN:-"Remove Mean Per voxel"}');
106. %data will be 60K x 200
107. 
      
108. %spacial ICA
109. [~,S,V]=svd(data, 'econ');
110. % V would be 200x200
111. 
      
112. Lambda=diag(S);
113. %Lambda = Lambda ./ sqrt(size(data, 1) - 1);
114. %[Lambda, inds] = sort(Lambda.^2); %make sv ascending
115. % no need to sort, default would be from Max to Min
116. 
      
117. LambdaSq  = Lambda.^2; %calculate variance for A'A
118. 
      
119. nc=${COMPNUM:-30}
120. evp=${EXPECTED_VAR_PERCENT:-0};
121. if evp>0,
122.     cumVar=100*cumsum(LambdaSq)/sum(LambdaSq);
123.     nc=find(cumVar>=evp, 1);
124. end
125. %by default, choose 30 components
126. pcaVar= sum(LambdaSq(1:nc))/sum(LambdaSq);
127. %but Variance Percent have higher priority
128. 
      
129. %dim reduction
130. Lambda=diag(Lambda(1:nc)); %vec to diag matrix
131. V=V(:,1:nc);
132. 
      
133. %map to PCA
134. whiteM = Lambda \ V'; %normalize by eigenvalue, complete whitening, will lose eigenvalue information
135. %30x200
136. 
      
137. dewhiteM = V * Lambda;
138. %200x30
139. 
      
140. pcasig=data*whiteM';  %60k x 30
141. pcasig=pcasig * diag(1./std(pcasig));
142. 
      
143. %print variance explained by the PCA
144. fprintf(2, '\n%d x %d, %d x %d ', size(data), size(pcasig));
145. fprintf(2, '\n### %.2g%% percent variance explained by %d PCA', 100*pcaVar, nc);
146. 
      
147. %MST
148. icasigR=cell(1, ${MST:-2});
149. for i=1:${MST:-2},
150.         [algo, W, A, icasigR{i}]=icatb_icaAlgorithm('${ICAALGORITHM:-ICA-EBM}',pcasig', ${OPT4ICA}, i);
151.         %input data should be 30x60k
152. end
153. [corrMetric, W, A, icasig, bestRun]=icatb_bestRunSelection(icasigR, pcasig');
154. %input data should be 30x60k
155. 
      
156. %back-proj loadings to be full, i.e., to make each subject have a loading
157. %pcasig = A*icasig (30x30 * 30x60k -> 30x60k)
158. 
      
159. loadings=dewhiteM * A; %200x30 * 30x30  -> 200x30
160. %data = loadings * icasig; %200x60k <- 200x30 * 30x60k
161. 
      
162. SM=icasig; tc=loadings;
163. icatb_saveICAData('$result', SM, tc, msk, size(SM, 1), hinfo, 'real');
164. dlmwrite('$resultTC', tc,'delimiter','\t')
165. 
      
166. % %rebuild 4D original dataset to nifti
167. % hinfo=spm_vol('$MASK');
168. % fname='${fpca}';
169. % for t=1:size(pcasig, 2),
170.     % odat=zeros(size(msko));
171.     % odat(msk)=pcasig(:, t);
172. 
      
173.     % hinfo.fname=fname;
174.     % hinfo.n(1)=t;
175.     % hinfo.dt(1)=16; % single / float32
176.     % spm_write_vol(hinfo, odat(:,:,:));
177. % end
178. catch
179.         xiaowei_lasterror;
180. end
181. debug=${DEBUG:-0};
182. if debug~=0,
183.     save('debug.mat');
184. end
185. 
      
186. exit
187. eom
188. 
      
189. matlab < $tmp
190. 
      
191. if IsFileLoadable $result; then
192. #cut Z-cat components back to each
193. if [[ 1 -lt $nnii ]]; then
194.     mkStack fix4sform
195.     #for((i=0;i < $nnii ;i++)); do
196.     for((i=$nnii-1;i>=0;i--)); do
197.         s=$(echo $i|awk -v az=$az '{print az*$1}')
198.         e=$(echo $i|awk -v az=$az '{print az*(1+$1)-1}')
199.         pops inputs nii #pop out in reverse order
200.         3dZcutup -prefix ${CUT_DIR:-.}/${nii}.nii -keep $s $e $result
201. 
      
202.         #save for nifti header sform/qform fixation, otherwise, the mm will be too huge for MNI
203.         pushs fix4sform ${CUT_DIR:-.}/${nii}.nii
204.     done
205.     #save nifti header of 1st dataset, which should have correct mm sform/qform matrix
206.     #then apply to remaining datasets
207.     pops fix4sform  nii  #save the 1st dataset's header
208.     elog "$nii sform/qform header will be saved"
209.     fslorient -getsform $nii > ${nii}.sform
210.     sform=${nii}.sform; onii=$nii
211.     for((i=$nnii-1;i>0;i--)); do
212.         pops fix4sform  nii #pop out in reverse order
213.         elog "fixing $nii header to the sform/qform of $onii"
214.         fslorient -setsform  $(cat $sform) $nii
215.         fslorient -copysform2qform $nii
216.     done
217. fi
218. else
219.     elog "Failed: $result "
220. fi
221. 
      
222. if [[ 0 -eq ${DEBUG:-0} ]]; then
223.     showStack ftmps
224.     eval "tmpf=$(showStack ftmps)"
225.     [[  -z "${tmpf[@]}"  ]] || rm ${tmpf[@]}
226.         rm $tmp
227. else
228.         elog $tmp
229. fi