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).

Input:
1a. (joint ICA over many components) a set of 3D NIFTI file
1b. (2nd level spatial ICA) a single 4D NIFTI file, while the temporal dimension made up from subjects' one corresponding spatial component across 2 groups, from back-reconstruction.
2. a mask in NIFTI format.
3. ICA dimension, i.e., how many components
4. MST number (default 10), i.e., how many runs of ICA to find a stable/reliable decomposition
Output:
1. a loading matrix (saved in csv and mat, each row corresponding to a subject)
2. a set of spatial maps (saved as NIFTI file), each component corresponding to a column of the loading matrix.
 
References:
A two-level ICA approach reveals important differences in the female brain response to thermal pain
X Song, S Bhinge, R Quiton, T Adali
ISBI 2018
 
Code:
$ cat ica.jointICA4nii.sh
  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. (Joint ICA assuming same loading matrix across joint components)
  22.     2. need manually concat in T dim if all inputs are 3D nifti files (Single spatial ICA but estimating across subjects for one spatial component)
  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
 
 
$ cat ica.jointICA4nii.sh
#!/usr/bin/env bash
#
# gicaS.sh
#% Init: 2013-11-06.13:37
#% Version: 20170614.17:28 Xiaowei Song
#% Copyright (C) 2013~2020 Xiaowei.Song <dawnwei.song@gmail.com>
#% http://restfmri.net
#% Distributed under terms of the AFL (Academy Free license).
#
CALLDIR="$( cd "$( dirname "${0}" )" && pwd )"
elog(){ echo "$@" 1>&2 ; }
 echo "$(date) | $0 $*" >> $(basename $0).LOG
source $(dawnbin)/slib.bash.sh
 
shopt -s expand_aliases; source $(dawnbin)/slib.stack.sh
usage(){ printf "
Usage: (%s)
    ${0} [options]
    options:
    Spacial ICA by using GIFT
    1. if multi nifti files input, default do Zconcat to make them in spatial dim
    2. need manually concat in T dim if all inputs are 3D nifti ( such as thickness ) files
" "$(grep -m1 '#% Version:' `readlink -f ${0}`)" 1>&2 ; }
while [ $# -gt 0 ]; do
    case "$1" in
        '-v') VERBOSE=1 ; shift ;;
        '-d') DEBUG=1 ; shift ;;
        '-m'|'--mask') MASK=$2 ; shift 2 ;;
        '-dm'|'--rm-mean') DEMEAN=$2; shift 2 ;;
        -k|'-n') COMPNUM=$2; shift 2 ;;
        #usually be 99->30
        -evp|--expected-var-percent) EXPECTED_VAR_PERCENT=$2; shift 2;;
        '-p'|'--prefix') PREFIX=$2; shift 2 ;;
                --mst) MST=$2; shift 2;;
        -rd|--cut-dir) CUT_DIR=$2; shift 2;;
        '-*') elog 'Unknown parameters'; usage; exit 0 ;;
        *) break ;;
    esac
done
verbose=${VERBOSE:-0};
 
if [ $# -eq 0 ] || ! IsFileExist $MASK ; then usage; exit 0; fi
nii=$1
bnii=$(basename $nii .nii)
#%parse ICA options
OPT4ICA=${ICAOPTIONS:-{}}
fpca="${PREFIX:-${bnii}}_pca.nii"
rm -f $fpca #clean, must do  since SPM8 will only overwrite certain time points
 
if ! IsFileLoadable $MASK; then
    elog "Mask not loadable: $MASK"
    #automask based on first input nii
    MASK=$(getmp).nii
    3dAutomask -prefix $MASK $1
fi
read ax ay az <<<$(3dinfo -ni -nj -nk $MASK)
 
 
mkStack inputs
mkStack ftmps
#build all input and zstack
lst=${PREFIX:-jointICA}.lst
while [[ $# -gt 0 ]]; do
    nii=$1
    pushs inputs $(rmExtDir4nii $nii)
    gz=${nii##*.}
    prefix=$(rmExtDir4nii $nii)
    if [[ gz == "$gz" ]]; then
        nii=$(getmp).nii
        pushs ftmps $nii
        gunzip -c $1 > $nii
    fi
    echo $(readlink -f $nii) ; shift;
done > $lst
showStack inputs
#exit
nnii=$(cat $lst|wc -l)
 
#build jICA input
nii=$(getmp).nii  ; pushs ftmps $nii
msk=$(getmp).nii ; pushs ftmps $msk
3dhdr $(cat $lst |xargs) $MASK
if [[ 1 -eq $nnii ]]; then
    3dcopy $(cat $lst |xargs) $nii #thus support nii.gz if only 1 input
    3dcopy $MASK $msk
else
    #3dZcat needs at least 2 inputs
    3dZcat -prefix $nii $(cat $lst |xargs)
    3dZcat -prefix $msk $(yes $MASK |head -n $nnii |xargs)
fi
 
 
result=${PREFIX:-${bnii}}_component_ica_.nii
resultTC=${PREFIX:-${bnii}}_timecourses_ica_.txt
rm $result $resultTC
 
tmp=$(getmp).m
cat > $tmp <<eom
try
 
[msko, hinfo]=icatb_read_data('$msk');
msk=find(msko);
 
data=icatb_read_data('$nii',[],msk);
data=icatb_preproc_data(data, '${DEMEAN:-"Remove Mean Per Timepoint"}');
%data will be 60K x 200
 
%spacial ICA
[~,S,V]=svd(data, 'econ');
% V would be 200x200
 
Lambda=diag(S);
%Lambda = Lambda ./ sqrt(size(data, 1) - 1);
%[Lambda, inds] = sort(Lambda.^2); %make sv ascending
% no need to sort, default would be from Max to Min
 
LambdaSq  = Lambda.^2; %calculate variance for A'A
 
nc=${COMPNUM:-30}
evp=${EXPECTED_VAR_PERCENT:-0};
if evp>0,
    cumVar=100*cumsum(LambdaSq)/sum(LambdaSq);
    nc=find(cumVar>=evp, 1);
end
%by default, choose 30 components
pcaVar= sum(LambdaSq(1:nc))/sum(LambdaSq);
%but Variance Percent have higher priority
 
%dim reduction
Lambda=diag(Lambda(1:nc)); %vec to diag matrix
V=V(:,1:nc);
 
%map to PCA
whiteM = Lambda \ V'; %normalize by eigenvalue, complete whitening, will lose eigenvalue information
%30x200
 
dewhiteM = V * Lambda;
%200x30
 
pcasig=data*whiteM';  %60k x 30
pcasig=pcasig * diag(1./std(pcasig));
 
%print variance explained by the PCA
fprintf(2, '\n%d x %d, %d x %d ', size(data), size(pcasig));
fprintf(2, '\n### %.2g%% percent variance explained by %d PCA', 100*pcaVar, nc);
 
%MST
icasigR=cell(1, ${MST:-2});
for i=1:${MST:-2},
        [algo, W, A, icasigR{i}]=icatb_icaAlgorithm('${ICAALGORITHM:-ICA-EBM}',pcasig', ${OPT4ICA}, i);
        %input data should be 30x60k
end
[corrMetric, W, A, icasig, bestRun]=icatb_bestRunSelection(icasigR, pcasig');
%input data should be 30x60k
 
%back-proj loadings to be full, i.e., to make each subject have a loading
%pcasig = A*icasig (30x30 * 30x60k -> 30x60k)
 
loadings=dewhiteM * A; %200x30 * 30x30  -> 200x30
%data = loadings * icasig; %200x60k <- 200x30 * 30x60k
 
SM=icasig; tc=loadings;
icatb_saveICAData('$result', SM, tc, msk, size(SM, 1), hinfo, 'real');
dlmwrite('$resultTC', tc,'delimiter','\t')
 
% %rebuild 4D original dataset to nifti
% hinfo=spm_vol('$MASK');
% fname='${fpca}';
% for t=1:size(pcasig, 2),
    % odat=zeros(size(msko));
    % odat(msk)=pcasig(:, t);
 
    % hinfo.fname=fname;
    % hinfo.n(1)=t;
    % hinfo.dt(1)=16; % single / float32
    % spm_write_vol(hinfo, odat(:,:,:));
% end
catch
        xiaowei_lasterror;
end
debug=${DEBUG:-0};
if debug~=0,
    save('debug.mat');
end
 
exit
eom
 
#runxvfb matlabg < $tmp
matlabj < $tmp
 
if IsFileLoadable $result; then
#cut Z-cat components back to each
if [[ 1 -lt $nnii ]]; then
    mkStack fix4sform
    #for((i=0;i < $nnii ;i++)); do
    for((i=$nnii-1;i>=0;i--)); do
        s=$(echo $i|awk -v az=$az '{print az*$1}')
        e=$(echo $i|awk -v az=$az '{print az*(1+$1)-1}')
        pops inputs nii #pop out in reverse order
        3dZcutup -prefix ${CUT_DIR:-.}/${nii}.nii -keep $s $e $result
 
        #save for nifti header sform/qform fixation, otherwise, the mm will be too huge for MNI
        pushs fix4sform ${CUT_DIR:-.}/${nii}.nii
    done
    #save nifti header of 1st dataset, which should have correct mm sform/qform matrix
    #then apply to remaining datasets
    pops fix4sform  nii  #save the 1st dataset's header
    elog "$nii sform/qform header will be saved"
    fslorient -getsform $nii > ${nii}.sform
    sform=${nii}.sform; onii=$nii
    for((i=$nnii-1;i>0;i--)); do
        pops fix4sform  nii #pop out in reverse order
        elog "fixing $nii header to the sform/qform of $onii"
        fslorient -setsform  $(cat $sform) $nii
        fslorient -copysform2qform $nii
    done
fi
else
    elog "Failed: $result "
fi
 
if [[ 0 -eq ${DEBUG:-0} ]]; then
    showStack ftmps
    eval "tmpf=$(showStack ftmps)"
    [[  -z "${tmpf[@]}"  ]] || rm ${tmpf[@]}
        rm $tmp
else
        elog $tmp
fi
AttachmentSize
slib.stack_.sh3.32 KB