FSL-MRS spec2nii twix conversion issue

jnlingad · May 27, 2022, 2:22am

Hi!

I’ve just begun using FSL-MRS (v1.1.12) and I’m running into some issues trying to convert my multiraid twix files acquired from a 3T Siemens Prisma with spec2nii (v0.4.2).

My multiraid file contains a “noise” dataset (-m 1):

spec2nii twix -v meas_MID01821_FID13842_Fus_Prep2_Ref2NoInline_SaveEach.dat -m 1
pymapVBVD version 0.4.8
Software version: VD
Contents of file meas_MID01821_FID13842_Fus_Prep2_Ref2NoInline_SaveEach.dat:  
Multiraid file, 2 files found.
Selecting file 1. Use -m option to change.
The file contains these evalinfo flags with dimensions and sizes as follows:
noise          :	Col, Cha, Lin, Ave  	[512  32 128   2]

and also my main “image” dataset (-m 2):

spec2nii twix -v meas_MID01821_FID13842_Fus_Prep2_Ref2NoInline_SaveEach.dat -m 2
pymapVBVD version 0.4.8
Software version: VD
Contents of file meas_MID01821_FID13842_Fus_Prep2_Ref2NoInline_SaveEach.dat:  
Multiraid file, 2 files found.
Selecting file 2. Use -m option to change.
The file contains these evalinfo flags with dimensions and sizes as follows:
image          :	Col, Cha, Ave, Phs  	[2080   32  150    2]

I was able to convert the main dataset “image”:

spec2nii twix -e image -m 2 -f Fus_run1 -j meas_MID01821_FID13842_Fus_Prep2_Ref2NoInline_SaveEach.dat
pymapVBVD version 0.4.8
Software version: VD
Converting twix file meas_MID01821_FID13842_Fus_Prep2_Ref2NoInline_SaveEach.dat.
Looking for evalinfo flag image.
Found data of size (2080, 32, 150, 2).                                       
Header extension validated!

However, I’m having trouble interpreting what it’s reading in here as the data size. The data are supposed to be 2048-point FIDs, collected on 32 channels, with 150 averages… but instead, the data size shown here shows that I have 2080-point FIDs, and I’m unclear about the 4th value of 2 represents here.

> pymapVBVD version 0.4.8
> Software version: VD
> 
> Converting twix file meas_MID01821_FID13842_Fus_Prep2_Ref2NoInline_SaveEach.dat.                                                                                                                   
> Looking for evalinfo flag image.
> Found data of size (2080, 32, 150, 2).                                                                                                                                                            
> Header extension validated!

Further, I’m unclear about the nature of the data stored in -m 1 has this datasize (512, 32, 128, 2), but I suspect this is my water unsuppressed data where 2 averages were acquired. I’d like to take a look at this after converting to .nii, but spec2nii fails:

spec2nii twix -e noise meas_MID01821_FID13842_Fus_Prep2_Ref2NoInline_SaveEach.dat -m 1 -f Fus_run1_wref
pymapVBVD version 0.4.8
Software version: VD
                                                                                                                                                                                                 Converting twix file meas_MID01821_FID13842_Fus_Prep2_Ref2NoInline_SaveEach.dat.                                                                                                                   
Looking for evalinfo flag noise.
Found data of size (512, 32, 128, 2).                                                                                                                                                             
Traceback (most recent call last):
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/bin/spec2nii", line 8, in <module>
    sys.exit(main())
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/spec2nii/spec2nii.py", line 509, in main
    spec2nii(*args)
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/spec2nii/spec2nii.py", line 260, in __init__
    args.func(args)
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/spec2nii/spec2nii.py", line 333, in twix
    self.imageOut, self.fileoutNames = process_twix(twixObj,
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/spec2nii/twixfunctions.py", line 95, in process_twix
    return process_svs(
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/spec2nii/twixfunctions.py", line 162, in process_svs
    meta_obj = extractTwixMetadata(twixObj['hdr'], basename(twixObj[dataKey].filename))
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/spec2nii/twixfunctions.py", line 440, in extractTwixMetadata
    return extractTwixMetadata_vx(mapVBVDHdr, original_file)
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/spec2nii/twixfunctions.py", line 567, in extractTwixMetadata_vx
    obj = nifti_mrs.hdr_ext(mapVBVDHdr['Meas'][('Frequency')] / 1E6,
KeyError: 'Frequency'

Lastly, I actually cannot view the main dataset that successfully converted to .nii successfully with mrs_tools:

mrs_tools vis --no_mean Fus_run1.nii.gz 
Performing coil combination
/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/fsl_mrs/utils/preproc/combine.py:46: RuntimeWarning:

divide by zero encountered in true_divide

/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/fsl_mrs/utils/preproc/combine.py:46: RuntimeWarning:

invalid value encountered in matmul

Traceback (most recent call last):
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/bin/mrs_tools", line 357, in <module>
    main()
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/bin/mrs_tools", line 129, in main
    args.func(args)
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/bin/mrs_tools", line 240, in vis
    vis_nifti_mrs(p)
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/bin/mrs_tools", line 177, in vis_nifti_mrs
    data = nifti_mrs_proc.coilcombine(data)
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/fsl_mrs/utils/preproc/nifti_mrs_proc.py", line 79, in coilcombine
    combinedc_obj[idx] = preproc.combine_FIDs(
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/fsl_mrs/utils/preproc/combine.py", line 158, in combine_FIDs
    return svd_reduce(FIDlist, W)
  File "/mnt/hippocampus/starkdata1/Jess/anaconda3/lib/python3.8/site-packages/fsl_mrs/utils/preproc/combine.py", line 69, in svd_reduce
    U, S, V = np.linalg.svd(FIDs, full_matrices=False)
  File "<__array_function__ internals>", line 5, in svd
  File "/home/jess/.local/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 1660, in svd
    u, s, vh = gufunc(a, signature=signature, extobj=extobj)
  File "/home/jess/.local/lib/python3.8/site-packages/numpy/linalg/linalg.py", line 97, in _raise_linalgerror_svd_nonconvergence
    raise LinAlgError("SVD did not converge")
numpy.linalg.LinAlgError: SVD did not converge

Has anyone run into this and might know how I could resolve these issues? Thanks so much in advance!

Jess

wclarke · June 2, 2022, 6:41am

Hi @jnlingad,

Thanks for this detailed bug report and for trying these tools. Sorry it’s taken me a few days to pick this up. I think there are a couple things going on here:

The 2080 vs 2048 time points is possibly expected. Sequences often add additional points to the ADC (read out event) before the centre of the echo is expected. This is so the artefacts seen at the beginning of the ADC don’t overlap the final data. Eventually you can discard these points with fsl_mrs_proc truncate ....
Regarding the two dimensions. Is the data either from a MEGA edited sequence, or is it a Siemens product sequence? Do you know what baseline you are using? If the latter, I think there’s a new set of product sequence which encodes a single water reference scan. This means that all the points, except a single average which contains the water reference in the first element of the final (size-2) AVE dimension are zeros. This is why the mrs_tools vis call fails.
I handle this for the new XA line of scanners (Vida etc) but perhaps this behaviour also occurs on some Prisma scanners.
Those noise scans in the first twix file of the multiraid -m 1 is effectively from another (non-spectroscopy) scanner calibration scan. That means spec2nii isn’t set up to handle them. However it might be worth looking into as the noise scans are useful.

jnlingad · June 11, 2022, 7:12am

Ah, I figured out that I can use mrs_tools split here to separate out the single water reference average from the water-suppressed data. Very excited to move forward with this. Thanks so much for your help!

Jess

greetingsproffalken · January 27, 2023, 12:49am

Running into this same issue on a Siemens. Can you elaborate on how you split? Thanks.

Helge · January 27, 2023, 2:09am

Hi @greetingsproffalken,

you can use mrs_tools split to get only the first average from the water reference file. For details see (https://open.win.ox.ac.uk/pages/fsl/fsl_mrs/data_handling.html) or type mrs_tools split -h in the terminal.

For example mrs_tools split --file water_ref.nii.gz --dim DIM_DYN --index 0 --output outputfolder generates two files one with only the first and one with the remaining entries along the DIM_DYN dimension.

Best,
Helge

jnlingad · January 27, 2023, 2:25am

Hi @greetingsproffalken,

Here’s my code for splitting MRS Siemens Twix (FSL-MRS version 2.0.2 was used here - not sure how this may have changed with the latest version):

        # split at the last dimension
        mrs_tools split \
            --file $OUTPUTPATH/${region}_${filenum}.nii.gz \
            --output $OUTPUTPATH \
            --dim DIM_USER_0 \
            --index 0 \
            --filename ${region}_${filenum}_split

        # the second file is the water-suppressed metabolite data
        mrs_tools info $OUTPUTPATH/${region}_${filenum}_split_2.nii.gz
        
        # the first is the water reference data: first we need to split that one at the 1 index
        mrs_tools split \
            --file $OUTPUTPATH/${region}_${filenum}_split_1.nii.gz \
            --output $OUTPUTPATH \
            --dim DIM_DYN \
            --index 1 \
            --filename $OUTPUTPATH/${region}_${filenum}_wref

        # check that the first file will be our water reference file which only consists of 2 transients
        mrs_tools info $OUTPUTPATH/${region}_${filenum}_wref_1.nii.gz
        
        # check that the second file contains the remaining transients of all zeros
        mrs_tools info $OUTPUTPATH/${region}_${filenum}_wref_2.nii.gz

        # delete split files that we don't need anymore
        rm $OUTPUTPATH/${region}_${filenum}_split_1.nii.gz
        rm $OUTPUTPATH/${region}_${filenum}_wref_2.nii.gz

You can then use mrs_tools vis on your outputs here with making sure to specify a broader ppm range with --ppmlim to inspect the water unsuppressed data.

greetingsproffalken · January 27, 2023, 2:28am

Appreciate the quick help from you both!

wclarke · January 30, 2023, 9:03am

Thanks @Helge and @jnlingad,

Just a note that mrs_tools has now been split off as a separate package called nifti-mrs. It still comes vendored with FSL-MRS but you can get it as standalone if you like to use it with other tools.

See GitHub - wtclarke/nifti_mrs_tools: Software tools for the NIfTI-MRS data format

greetingsproffalken · February 10, 2023, 10:51pm

Thanks again all. So I finally had time to loop in on this. So, I first generated the nifti files from my raw dat files (2 scans: water suppressed and water ref) with spec2nii. I used -m2 option to convert the non-noise files. The water supp scan looked as follows:

[‘DIM_COIL’, ‘DIM_DYN’, ‘DIM_USER_O’]
(1,1,1,2080,32,128,2)

So I then ran mrs_tools split on that file with options --dim DIM_USER_0 --index 0 which did remove the ‘2’ dim. But, it gave me files tagged high and low as outputs. I am using mrs_tools version 2.0.9. The high tagged file is a much larger file, but they look the same with mrs_tools info, and I can only run mrs_tools vis on the high file, and I get a reasonable looking output.

What do the high and low files correspond to? Do I still need to split 2 averages out from the high file (i.e. using --dim DIM_DYN --index 1), or could I just use the high file and the water ref scan that we ran (high file, split as above) to proceed with preprocessing like I learned at the Oxford FSL MRS course? Lastly, do you recommend installing the nifti-mrs package moving forward? Thanks, I’m a rook!

wclarke · February 11, 2023, 7:15pm

It sounds like the ‘low’ file contains just (or elements with only) zeros. That would explain why it is a) smaller, and b) I presume mrs_tools vis fails on the coil combination step (the SVD fails on all zero data, something I need to add nicer error handling of). I.e. of the 128 remaining dynamic elements there some are empty. If you have a suspicion of which ones you can run mrs_tools split again on the other dimension (DIM_DYN) to split things down further to remove the zeros.

If you have a new version of FSL-MRS you will have the package nifti-mrs package (they come vendored together). It might be easiest to load the file in to a python environment (e.g. a notebook) and use the nifti-mrs package to have a play with the data and see what is in there.

wclarke · February 11, 2023, 7:18pm

Clearly there needs to be a) work from sequence developers to not be lazy in their data labelling (the EVALINFO is there to stop these issues), and b) I need to include a tool in mrs_tools that can deal with empty transients.