MRSI ROI selection

This will ideally provide a graphical interface to guide placing coordinates and extracting MR spectrum – a python port of the slow-to-load, hard-to-hack, and unusable-over-X11Forwarding matlab version in matlab/.

Citation

This and QPASA slice warping used first in Perica et al, 2022

Regions of interest (ROIs) were defined in MNI-space using a custom ROI atlas, and included right and left dorsolateral prefrontal cortex (DLPFC), right and left anterior insula (AI), anterior cingulate cortex (ACC), and medial PFC (MPFC) (Fig. 1b). An automated nonlinear registration using FSL’s FNIRT was performed to map those ROIs into each subject’s native space. ROI centers were visually inspected and manually adjusted for each subject with in-house developed software in order to maximize spectral quality and gray matter content of the voxel. To ensure manual adjustment did not move ROI centers outside of the region of interest, the adjusted ROIs were registered back to standard space where their position was confirmed with the Eickhoff-Zielles maximum probability map of the Talairach-Tournoux atlas using AFNI’s whereami (Cox, 1996). The ROI coordinates were then used in the anatomically guided voxel-shifting reconstruction of MRSI data (Hetherington et al., 2007) to gather spectroscopic data over the ROI

Usage

## SETUP. clone and get depends. only need to run once
# get MRSIcoord.py code
git clone https://github.com/LabNeuroCogDevel/MRSIcoord.py && cd MRSIcoord.py
# get any missing python libraries
pip install -r requirements.txt

## run GUI with example input data
./grid.py \
   -s test/data/siarray.1.1  -r test/data/rorig.nii --gm_mask test/data/gm_sum.nii.gz \
   -i test/data/roi_pos.txt  --rois roi4 roi6 

Inputs

-s siarray.1.1	spectroscipy data. 1024 complex (1024 real, 1024 imaginary) points.
-r rorig.nii	structural MR registered to scout position
--gm_mask gm_sum.nii	binary mask for if voxel is gm. sum of id Freesurfer gray matter regions, also coregistered to the scout
-i coords.txt	initial guess at ROI positions (from mni atlas warp, center of mass). tab delim w/columns: roilabel, x, and y

Output

spectrum.xx.yy	reconstructed spectrum (2D FFT) at coordinate xx, yy. One file for each ROI.
spectrum.xx.yy.dir/spreedsheet.csv	lcmodel’ed metabolite concentrations

Pipeline

• read inputs
• place each ROI (get coordinates)
  • optimize for spectrum quality (e.g. far enough from skull)
  • greatest GM count
  • still classified as region of interest (anatomy)
• generate spectrum files
• run lcmodel. (to create binary exe, see make lcmodel/lcmodel)

Notes

• Defaults assume an 216x216 1mm T1w/anatomical resolution and 9x9x10mm MRSI slice
• basis function specific to MRRC 7T (PFC and Hc slices) ./lcmodel/basis-sets/gamma_7TJref_te34_297mhz_no2HG_1.basis
• need to compile lcmodel (https://github.com/schorschinho/LCModel) for your system and put the binary in the lcmodel folder as lcmodel/lcmodel. there is a Makefile rule for this
• “sid3” coordinate/matrix order convention is inherited from MRRC used for kspace+lcmodel: COL, ROW? (fortran column-major)
  • --roi_initial roi.txt is like (label, x, y); sid is like (y, 216-x)
  • output directories like specturm.yy.xx are labeled with sid3 coordinates
• siarray.1.1 is from VB/VD MRI raw Siemens ‘twix’ .dat. See pyMapVBVD or original matlab. Not created here.

Standalone python scripts

see Makefile for examples in use

spectrum

mkspectrum extracts single series (1024 complex values) around a given coordinate (from 216x216 grid)

./mkspectrum test/data/siarray.1.1 216 --pos 112 88 out/
# creates out/spectrum.112.88

lcmodel

lcmodel.py can be used as a stand alone script to run lcmodel on the above

./lcmodel.py out/spectrum.112.88
# creates  test/data/spectrum.112.88.dir/spreadsheet.csv

Porting From Matlab

SID3 coordinate/matrix order

TODO:

• python saves coord and sid_coord files (need to replace LNCD pipeline)

# converting sid3 spectrum.xx.yy filename to match afni's center of mass warped rois '3dCM -local_ijk'
sid3(){ awk '{print 216-$3+1 "\t" 216-$2+1}' "$@"; }

roi_slice_ijk(){
  # incomplete summary of /Volumes/Hera/Projects/7TBrainMech/scripts/mri/MRSI_roi/000_setupdirs.bash
  # (1) warp mni atlas roi to slice space. (2) keep only center slice (match MRSI acq). (3) get center of rois
  # roi centers used as starting point for gui placement
  applywarp -o $outimg -i $mni_atlas -r $pfc_ref -w $mni_to_t1 --postmat=$t1_to_pfc --interp=nn
  3dcalc -a "$outimg" -expr "equals(k,$slice_num_0)*a" -prefix middle_slice.nii.gz -overwrite
  3dresample -prefix $res_img -inset middle_slice.nii.gz -master rorig.nii.gz
  3dCM -local_ijk -all_rois $res_img | egrep '^[0-9]|#ROI'|paste - - |cut -f2-4 -d" "
}

Read/Write isn’t lossless?

After IFFT, matlab code saves to kspace.1.1. Reading this file back in and comparing to itself we can see lossy-ness around 10^-3. The figure max color is 10^-5. Plot shows siarray (python) IFFT against itself (matlab fwrite exported version). The same difference is seen comparing fwrite output with matlab’s IFFT.

   import matplotlib.pyplot as plt
   from siarray import SIArray
   # calc data
   SI = SIArray('test/data/siarray.1.1')
   SI.IFFTData()
   # read stored (matlab fwrite)
   with open('test/data/matlab/kspace.1.1', 'r') as f:
	 kspace = np.fromfile(f, '<4f')

   reread = kspace.reshape(24**2, 1024*2).T 
   orig = SI.kspace.reshape(24**2, 1024*2).T
   # see difference
   plt.imshow(abs(orig - reread))
   plt.clim([0,10**-5])
   plt.savefig('imgs/lossy-kspace.png', bbox_inches='tight')

Comparing python and matlab

test/genrate_mat.m runs through the spectrum pipeline and saves out matfiles to test the python code against.

python -m pytest

Example data from

find /Volumes/Hera/Projects/7TBrainMech/subjs/11743_20190802/slice_PFC/MRSI_roi/ -maxdepth 2 -iname 'rorig.nii' -or -iname 'mprage_middle.mat' -or -iname 'siarray.1.1' |
 xargs -I{} cp {} test/data/

Testing MATLAB

spectrum lossy read/write

Same gen_spectrum is off by at most .0023 on a value of 5422.2

si='/Volumes/Hera/Projects/7TBrainMech/subjs/10129_20180917/slice_PFC/MRSI_roi/raw/siarray.1.1'
gen_spectrum(si, 216, [112, 104], '/tmp')
a=fread(fopen('/Volumes/Hera/Projects/7TBrainMech/subjs/10129_20180917/slice_PFC/MRSI_roi/raw/spectrum.112.104'),'float')
b=fread(fopen('/tmp/spectrum_112.104'),'float');

% exactly correlated
corr(a,b) % 1.0000

% but not identical (off by .0229 on value of ~5000)
[v,i] = max(abs(a-b)); v, a(i), b(i),
%   0.0229
%   5.4422e+03
%   5.4421e+03

Positions

% ./grid.py  -s test/data/siarray.1.1  -r test/data/rorig.nii  -i test/data/roi_pos.txt  --rois roi4 roi6
cd matlab
coord_mover('test', 'subjcoords', '../test/data/pos_z.txt', 'brain', '../test/data/rorig.nii')
fid = fopen('../test/data/WF/spectrum.78.66'); ml7866 = fread(fid,'single');
fid = fopen('../out/spectrum.78.66'); py7866 = fread(fid,'single');
hist(ml7866 - py7866)
max(abs(ml7866 - py7866))
%   9.3126e-04
mean(abs(ml7866 - py7866))
%   4.0350e-05
[std(py7866), max(py7866)]
%   1.0e+03 * 0.2648    1.6928