Skip to content

GitLab

Commit ca2b1170 authored by William Clarke's avatar William Clarke
Browse files

Tidying up. Add quantificcation info to report.

parent 88d872a5
Hide whitespace changes
Inline Side-by-side
README.md
View file @ ca2b1170
......@@ -8,7 +8,9 @@ FSL-MRS is a collection of python modules and wrapper scripts for pre-processing
---
### Installation
git clone https://git.fmrib.ox.ac.uk/saad/fsl_mrs.git
To get the packaged example data, make sure [git-lfs](https://git-lfs.github.com/) is installed.
git clone --recurse-submodules https://git.fmrib.ox.ac.uk/saad/fsl_mrs.git
cd fsl_mrs
pip install .
......
docs/user_docs/install.rst
View file @ ca2b1170
......@@ -5,9 +5,13 @@ FSL-MRS can currently be installed using two methods
1. From GitLab
~~~~~~~~~~~~~~
Download or clone from |fsleyes_gitlab|_.
Download or clone from |fsleyes_gitlab|_. `Git LFS <https://git-lfs.github.com/>`_ must be installed to download package data.
Run `pip install .` in the top level directory once downloaded.
::
git clone --recurse-submodules https://git.fmrib.ox.ac.uk/saad/fsl_mrs.git
cd fsl_mrs
pip install .
2. From Conda
~~~~~~~~~~~~~
......
fsl_mrs/core/MRS.py
View file @ ca2b1170
......@@ -400,8 +400,9 @@ class MRS(object):
Parameters
----------
ppmlist : default is [0.9,1.7,1.4,1.2,2.0]
ppmlist : default is [0.9,1.2,1.4,1.7,[2.08,2.25,1.95,3.0]]
amplist : default is [3.0,2.0,2.0,2.0,[1.33,0.33,0.33,0.4]]
gamma,sigma : float parameters of Voigt blurring
"""
if ppmlist is None:
......
fsl_mrs/scripts/fsl_mrs
View file @ ca2b1170
......@@ -295,8 +295,6 @@ def main():
# Echo time
if args.TE is not None:
echotime = args.TE*1E-3
elif 'TE' in dataheader:
echotime = dataheader['TE']
elif 'meta' in basisheader:
if 'TE' in basisheader['meta']:
echotime = basisheader['meta']['TE']
......@@ -304,6 +302,8 @@ def main():
echotime *= 1E-3
else:
echotime = None
elif 'TE' in dataheader:
echotime = dataheader['TE']
else:
echotime = None
# Internal and Water quantification if requested
......@@ -339,23 +339,20 @@ def main():
print('--->> Saving output files to {}\n'.format(args.output))
res.to_file(filename=os.path.join(args.output,'results_table.csv'),what='concentrations')
res.to_file(filename=os.path.join(args.output,'summary.csv'),what='summary')
res.to_file(filename=os.path.join(args.output,'concentrations.csv'),what='concentrations')
res.to_file(filename=os.path.join(args.output,'qc.csv'),what='qc')
res.to_file(filename=os.path.join(args.output,'all_parameters.csv'),what='parameters')
# Save image of MRS voxel
location_fig = None
if args.t1 is not None:
datatype = mrs_io.check_datatype(args.data)
if datatype == 'NIFTI':
fig = plotting.plot_world_orient(args.t1,args.data)
fig.tight_layout()
fig.savefig(os.path.join(args.output,'voxel_location.png'),
bbox_inches='tight',facecolor='k')
# Create short HTML report
#fig = plotting.plotly_fit(mrs,res,ppmlim=ppmlim,proj='abs')
#plotly.io.write_html(fig, file=os.path.join(args.output,'short_report.html'))
location_fig = os.path.join(args.output,'voxel_location.png')
fig.savefig(location_fig,bbox_inches='tight',facecolor='k')
# Save quick summary figure
report.fitting_summary_fig(mrs,res,
......@@ -368,8 +365,9 @@ def main():
filename=os.path.join(args.output,'report.html'),
fidfile=args.data,
basisfile=args.basis,
h2ofile=args.h2o,
date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M"))
h2ofile=args.h2o,
date=datetime.datetime.now().strftime("%Y-%m-%d %H:%M"),
location_fig = location_fig)
......
fsl_mrs/utils/quantify.py
View file @ ca2b1170
......@@ -46,7 +46,7 @@ def quantifyWater(mrs,H2OFID,refFID,referenceName,concentrations,metaboliteNames
Args:
mrs (MRS obj): Current MRS object
H2OFID (np.array): FID of wate reference
H2OFID (np.array): FID of water reference
refFID (np.array): FID of fitted reference metabolite
referenceName (str): Name of reference metabolite
concentrations (np.array): All metabolite raw concentrations
......@@ -57,15 +57,16 @@ def quantifyWater(mrs,H2OFID,refFID,referenceName,concentrations,metaboliteNames
verbose (bool): Verbose output
Returns:
conc_molal (float): Scaling parameter to convert raw fitted concnetrations to molality units of mols/kg
conc_molar (float): Scaling parameter to convert raw fitted concnetrations to molarity units of mols/dm^3
conc_molal (float): Scaling parameter to convert raw fitted concentrations to molality units of mols/kg
conc_molar (float): Scaling parameter to convert raw fitted concentrations to molarity units of mols/dm^3
dict : Quantification information
"""
# Calculate observed areas
SMObs = calculate_area(mrs,refFID,ppmlim=reflimits)
SH2OObs = calculate_area(mrs,H2OFID,ppmlim=None)
# Calculate concnetration scalings
# Calculate concentration scalings
# EQ 4 and 6 in https://doi.org/10.1002/nbm.4257
# conc_molal = (SMObs *(Q.f_GM*Q.R_H2O_GM + Q.f_WM*Q.R_H2O_WM + Q.f_CSF*Q.R_H2O_CSF)/(SH2OObs*(1-Q.f_CSF)*Q.R_M)) \
# * (H2O_PROTONS/refProtons)\
......@@ -98,7 +99,7 @@ def quantifyWater(mrs,H2OFID,refFID,referenceName,concentrations,metaboliteNames
print(f'Final molality scaling = {conc_molal:0.2e}')
print(f'Final molarity scaling = {conc_molar:0.2e}')
return conc_molal,conc_molar
return conc_molal,conc_molar,{'q_info':Q, 'metab_amp':SMObs, 'water_amp':SH2OObs}
class QuantificationInfo(object):
""" Class encapsulating the information required to run internal water quantification scaling."""
......@@ -111,7 +112,9 @@ class QuantificationInfo(object):
tissueFractions (dict:float): Tissue volume fractions, must contain 'WM', 'GM', 'CSF' fields.
tissueDensity (dict:float, optional): Tissue volume fractions, must contain 'WM', 'GM', 'CSF' fields. In units of g/ml.
"""
self.TE = TE
self.T2 = T2
if tissueFractions is not None:
if tissueDensity is None:
self.d_GM = TISSUE_WATER_DENSITY['GM']
......
fsl_mrs/utils/report.py
View file @ ca2b1170
......@@ -19,6 +19,7 @@ from fsl_mrs.utils import misc
import nibabel as nib
import base64
# class Report(object):
......@@ -133,6 +134,21 @@ def create_plotly_div(mrs,res):
fig = plotting.plot_indiv_stacked(mrs,res,ppmlim=res.ppmlim)
divs['indiv'] = to_div(fig) #plotly.offline.plot(fig, output_type='div',include_plotlyjs='cdn')
# Quantification table
if (res.concScalings['molality'] is not None) and hasattr(res.concScalings['info']['q_info'],'f_GM'):
Q = res.concScalings['info']['q_info']
quant_df = pd.DataFrame()
quant_df['Tissue-water densities'] = [Q.d_GM,Q.d_WM,Q.d_CSF]
quant_df['Tissue volume fractions'] = [Q.f_GM,Q.f_WM,Q.f_CSF]
quant_df['Water T2s'] = [Q.T2['H2O_GM'],Q.T2['H2O_WM'],Q.T2['H2O_CSF']]
quant_df.index = ['GM', 'WM', 'CSF']
quant_df.index.name = 'Tissue'
quant_df.reset_index(inplace=True)
tab = plotting.create_table(quant_df)
fig = go.Figure(data=[tab])
fig.update_layout(height=100,margin=dict(l=0,r=0,t=0,b=0))
divs['quant_table'] = to_div(fig)
return divs
......@@ -239,7 +255,7 @@ def static_image(imgfile):
# fig.save(os.path.join(outdir,'voxplot.png'))
# return
def create_report(mrs,res,filename,fidfile,basisfile,h2ofile,date):
def create_report(mrs,res,filename,fidfile,basisfile,h2ofile,date,location_fig = None):
divs= create_plotly_div(mrs,res)
......@@ -283,7 +299,8 @@ def create_report(mrs,res,filename,fidfile,basisfile,h2ofile,date):
<a href="#qc">QC</a> -
<a href="#uncertainty">Uncertainty</a> -
<a href="#realimag">Real/Imag</a> -
<a href="#metabs">Metabs</a> -
<a href="#metabs">Metabs</a> -
<a href="#quantification">Quantification</a> -
<a href="#methods">Methods</a>
</center>
<hr>
......@@ -299,12 +316,9 @@ def create_report(mrs,res,filename,fidfile,basisfile,h2ofile,date):
template+=section
# Location
voxpngfile=os.path.join(os.path.split(filename)[0],"voxel_location.png")
#voxplothtml=f'<h1>Voxel location</h1><img src="{voxpngfile}" alt="Voxel location" width="700"></img><hr>'
#template+=voxplothtml
fig = static_image(voxpngfile)
template+=f"<h1>Voxel location</h1><div>{to_div(fig)}</div><hr>"
if location_fig is not None:
fig = static_image(location_fig)
template+=f"<h1>Voxel location</h1><div>{to_div(fig)}</div><hr>"
# Tables section
section=f"""
......@@ -347,12 +361,49 @@ def create_report(mrs,res,filename,fidfile,basisfile,h2ofile,date):
"""
template+=section
# Quantification information
# Table of CSF,GM,WM
# Fractions
# T2 info (water)
# T2 info (metab)
# Tissue water densities
# Relaxation corrected water
# Relaxation correction for metab
# Final scalings for molality & molarity
if res.concScalings['molality'] is not None:
Q = res.concScalings['info']['q_info']
relax_water_conc = Q.relaxationCorrWater_molar
metabRelaxCorr = Q.relaxationCorrMetab
if hasattr(res.concScalings['info']['q_info'],'f_GM'):
section=f"""
<h1><a name="quantification">Quantification information</a></h1>
{divs['quant_table']}
<p>The metabolite T<sub>2</sub> is {1000*Q.T2['METAB']} ms.<br>
The sequence echo time (T<sub>E</sub>) is {1000*Q.TE} ms.<br>
<p>The T<sub>2</sub> relaxation corrected water concentration is {relax_water_conc:0.0f} mmol/kg.<br>
The metabolite relaxation correction (1/e<sup>(-T<sub>E</sub>/T<sub>2</sub>)</sup>) is {metabRelaxCorr:0.2f}.<br>
<p>The final raw concentration to molarity and molality scalings are {res.concScalings["molarity"]:0.2f} and {res.concScalings["molality"]:0.2f}.
<hr>
"""
else:
section=f"""
<h1><a name="quantification">Quantification information</a></h1>
<p>The water T<sub>2</sub> is {1000*Q.T2['H2O']} ms.<br>
The metabolite T<sub>2</sub> is {1000*Q.T2['METAB']} ms.<br>
The sequence echo time (T<sub>E</sub>) is {1000*Q.TE} ms.<br>
<p>The T<sub>2</sub> relaxation corrected water concentration is {relax_water_conc:0.0f} mmol/kg.<br>
The metabolite relaxation correction (1/e<sup>(-T<sub>E</sub>/T<sub>2</sub>)</sup>) is {metabRelaxCorr:0.2f}.<br>
<p>The final molarity and molality scalings are {res.concScalings["molarity"]:0.2f} and {res.concScalings["molality"]:0.2f}.
<hr>
"""
template+=section
# Details of the methods
#methodsfile="/path/to/file"
#methods = np.readtxt(methodsfile)
from fsl_mrs import __version__
methods=f"<p>Fitting of the SVS data was performed using a Linear Combination model as described in [1] and as implemented in FSL-MRS version {__version__}, part of FSL (FMRIB's Software Library, www.fmrib.ox.ac.uk/fsl). Briefly, basis spectra are fitted to the complex FID in the frequency domain. The basis spectra are shifted and broadened with parameters fitted to the data and grouped into {max(res.metab_groups)+1} metabolite groups. A complex polynomial baseline is also concurrrently fitted (order={res.baseline_order}). Model fitting was performed using the {res.method} algorithm.<p><h3>References</h3><p>[1] Clarke W, Jbabdi S. FSL-MRS: An end-to-end spectroscopy analysis package (2020)."
methods=f"<p>Fitting of the SVS data was performed using a Linear Combination model as described in [1] and as implemented in FSL-MRS version {__version__}, part of FSL (FMRIB's Software Library, www.fmrib.ox.ac.uk/fsl). Briefly, basis spectra are fitted to the complex FID in the frequency domain. The basis spectra are shifted and broadened with parameters fitted to the data and grouped into {max(res.metab_groups)+1} metabolite groups. A complex polynomial baseline is also concurrrently fitted (order={res.baseline_order}). Model fitting was performed using the {res.method} algorithm.<p><h3>References</h3><p>[1] Clarke WT, Jbabdi S. FSL-MRS: An end-to-end spectroscopy analysis package (2020)."
section=f"""
<h1><a name="methods">Analysis methods</a></h1>
<div>{methods}</div>
......
fsl_mrs/utils/results.py
View file @ ca2b1170
......@@ -29,7 +29,7 @@ class FitRes(object):
self.baseline_order = baseline_order
self.base_poly = B
self.concScalings={'internal':None,'internalRef':None,'molarity':None,'molality':None}
self.concScalings={'internal':None,'internalRef':None,'molarity':None,'molality':None,'info':None}
def loadResults(self,mrs,fitResults):
"Load fitting results and calculate some metrics"
......@@ -111,7 +111,7 @@ class FitRes(object):
Q = quant.loadDefaultQuantificationInfo(TE,tissueFractions,mrs.centralFrequency/1E6)
else:
Q = quant.QuantificationInfo(TE,T2,tissueFractions)
molalityScaling,molarityScaling = quant.quantifyWater(mrs,
molalityScaling,molarityScaling,quant_info = quant.quantifyWater(mrs,
waterRefFID,
refFID,
waterReferenceMetab,
......@@ -122,9 +122,9 @@ class FitRes(object):
reflimits=reflimits,
verbose=verbose)
self.concScalings={'internal':internalRefScaling,'internalRef':self.intrefstr,'molarity':molarityScaling,'molality':molalityScaling}
self.concScalings={'internal':internalRefScaling,'internalRef':self.intrefstr,'molarity':molarityScaling,'molality':molalityScaling,'info':quant_info}
else:
self.concScalings={'internal':internalRefScaling,'internalRef':self.intrefstr,'molarity':None,'molality':None}
self.concScalings={'internal':internalRefScaling,'internalRef':self.intrefstr,'molarity':None,'molality':None,'info':None}
def combine(self,combineList):
"""Combine two or more basis into single result"""
......@@ -247,35 +247,85 @@ class FitRes(object):
def to_file(self,filename,what='concentrations'):
"""
Save results to a csv file
Save results to csv file
Parameters:
-----------
filename : str
what : one of 'concentrations, 'qc', 'parameters'
what : one of 'summary', 'concentrations, 'qc', 'parameters'
"""
df = pd.DataFrame()
if what == 'concentrations':
if what == 'summary':
df = pd.DataFrame()
df['Metab'] = self.metabs
df['Raw conc'] = self.getConc()
if self.concScalings['internal'] is not None:
concstr = f'/{self.concScalings["internalRef"]}'
df[concstr] = self.getConc(scaling='internal')
else:
df['Raw conc'] = self.getConc()
if self.concScalings['molality'] is not None:
df['mMol/kg'] = self.getConc(scaling='molality')
if self.concScalings['molarity'] is not None:
df['mM'] = self.getConc(scaling='molarity')
df['%CRLB'] = self.perc_SD[:self.numMetabs]
df['%CRLB'] = self.getUncertainties()
SNR = self.getQCParams()[0]
SNR.index = SNR.index.str.replace('SNR_','')
SNR.index.name = 'Metab'
SNR = SNR.reset_index(name='SNR')
df = df.merge(SNR,how='outer')
FWHM = self.getQCParams()[1]
FWHM.index = FWHM.index.str.replace('fwhm_','')
FWHM.index.name = 'Metab'
FWHM = FWHM.reset_index(name='FWHM')
df = df.merge(FWHM,how='outer')
elif what == 'concentrations':
scaling_type = ['raw']
if self.concScalings['internal'] is not None:
scaling_type.append('internal')
if self.concScalings['molality'] is not None:
scaling_type.append('molality')
if self.concScalings['molarity'] is not None:
scaling_type.append('molarity')
all_df = []
for st in scaling_type:
df = self.getConc(scaling=st,function=None).T
df.insert(0,'mean',df.mean(axis=1))
df.insert(1,'standard deviation',self.getUncertainties(type=st))
df.insert(0,'scaling',st)
df.index.name = 'metabolite'
all_df.append(df)
df = pd.concat(all_df)
df.reset_index(inplace=True)
elif what == 'qc':
SNR = self.SNR.peaks.T
SNR.index = SNR.index.str.replace('SNR_','')
SNR.insert(0,'mean',SNR.mean(axis=1))
SNR.insert(0,'Parameter','SNR')
SNR
FWHM = self.FWHM.T
FWHM.index = FWHM.index.str.replace('fwhm_','')
FWHM.insert(0,'mean',FWHM.mean(axis=1))
FWHM.insert(0,'Parameter','FWHM')
df = pd.concat([SNR,FWHM])
df.index.name = 'metabolite'
df.reset_index(inplace=True)
elif what == 'qc':
pass
# snr,fwhm = self.getQCParams()
# df['Measure'] = ['SNR']
# df['Value'] = [self.snr]
elif what == 'parameters':
df['Name'] = self.fitResults.columns
df['Value'] = self.fitResults.mean().to_numpy()
df = self.fitResults.T
df.insert(0,'mean',self.fitResults.T.mean(axis=1))
df.insert(1,'std',self.fitResults.T.std(axis=1))
df.index.name = 'parameter'
df.reset_index(inplace=True)
df.to_csv(filename,index=False,header=True)
......@@ -332,7 +382,7 @@ class FitRes(object):
if phi0.lower() == 'degrees':
p0 *= np.pi/180.0
elif phi0.lower() == 'radians':
elif (phi0.lower() == 'radians') or (phi0.lower() == 'raw'):
pass
else:
raise ValueError('phi0 must be degrees or radians')
......@@ -343,8 +393,10 @@ class FitRes(object):
p1 *= -180.0/np.pi * self.hzperppm
elif phi1.lower() == 'deg_per_hz':
p1 *= -180.0/np.pi * 1.0
elif phi1.lower() == 'raw':
pass
else:
raise ValueError('phi1 must be seconds or deg_per_ppm or deg_per_hz ')
raise ValueError('phi1 must be seconds, deg_per_ppm, deg_per_hz or raw.')
return p0,p1
......@@ -364,8 +416,10 @@ class FitRes(object):
shiftParams *= 1/(2*np.pi*self.hzperppm)
elif units.lower() == 'hz':
shiftParams *= 1/(2*np.pi)
elif units.lower() == 'raw':
pass
else:
raise ValueError('Units must be Hz or ppm.')
raise ValueError('Units must be Hz, ppm or raw.')
return shiftParams
......@@ -386,8 +440,10 @@ class FitRes(object):
gamma *= 1/(np.pi)
elif units.lower() == 'ppm':
gamma *= 1/(np.pi*self.hzperppm)
elif units.lower() == 'raw':
pass
else:
raise ValueError('Units must be Hz or ppm.')
raise ValueError('Units must be Hz, ppm or raw.')
combined = gamma
return combined,gamma
elif self.model=='voigt':
......@@ -412,8 +468,10 @@ class FitRes(object):
elif units.lower() == 'ppm':
gamma *= 1/(np.pi*self.hzperppm)
sigma *= 1/(np.pi*self.hzperppm)
elif units.lower() == 'raw':
pass
else:
raise ValueError('Units must be Hz or ppm.')
raise ValueError('Units must be Hz, ppm or raw.')
combined = gamma/2 + np.sqrt((gamma**2/4.0)+(sigma*2*np.log(2))**2)
return combined,gamma,sigma
......@@ -445,7 +503,7 @@ class FitRes(object):
def getUncertainties(self,type='percentage',metab=None):
""" Return the uncertainties on concentrations.
""" Return the uncertainties (SD) on concentrations.
Can either be in raw, molarity or molality or percentage uncertainties.
"""
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment