{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## `T1`, `MTR` and `MTsat` in-vivo\n",
"\n",
"\n",
"\n",
"3 healthy participants were scanned with the following `MTsat` protocol: \n",
"\n",
"```{admonition} Click here to reveal the MTsat protocol.\n",
":class: dropdown\n",
"\n",
"|Parameter (PDw/MTw/T1w)|G1NATIVE|S1NATIVE|S2NATIVE| VENUS|\n",
"|---------|-----------|-----------|-----------|-----|\n",
"| Sequence Name | SPGR | FLASH| FLASH| mt_sat v1.0 |\n",
"| Flip Angle (°) | 6/6/20 | 6/6/20| 6/6/20| 6/6/20 |\n",
"| TR (ms) | 32/32/18 | 32/32/18| 32/32/18| 32/32/18 |\n",
"| TE (ms) | 4 | 4| 4| 4|\n",
"| FOV (cm) | 25.6 | 25.6| 25.6| 25.6|\n",
"| Acquisition Matrix | 256x256 | 256x256| 256x256| 256x256|\n",
"| Receiver BW (kHz) | 62.5 | 62.5| 62.5| 62.5|\n",
"| RF Phase Increment (°) | 115.4 | 50| 50| 117|\n",
"| MT Pulse | off/on/off | off/on/off| off/on/off| off/on/off |\n",
"| MT Pulse Duration (ms) | 8 | 10| 10| 12 |\n",
"| MT Pulse Shape | Fermi | Gaussian| Gaussian| Fermi |\n",
"| MT Pulse Offset (Hz) | 1200 | 1200| 1200| 1200 |\n",
"\n",
"```\n",
"\n",
"The following code cell loads all the necessary Python packages, defines common functions and downloads the data to `/tmp/VENUS` directory. If data already exists, the download will be skipped.\n",
"\n",
"The number of data points per kernel-density estimations (KDE) were limited to 10k for faster visualization. Lower and upper limits of the distributions are set according to the original publication. \n",
"\n",
"\n",
"`````{admonition} Tip\n",
":class: tip\n",
"* You can use the legend on the right to display a selection of KDEs. For example, double click on G1VENUS to single out vendor-neutral distributions acquired on the GE scanner. \n",
"* At the top right corner of the plot background, you can switch the hover mode to `compare data on hover` to compare KDEs point by point.\n",
"`````"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"import plotly.graph_objs as go\n",
"import plotly.express as px\n",
"from plotly import tools, subplots\n",
"import numpy as np\n",
"import pandas as pd\n",
"import ipywidgets as widgets\n",
"import math\n",
"from plotly.offline import download_plotlyjs, init_notebook_mode, plot, iplot\n",
"from ipywidgets import interact, interactive, fixed, interact_manual\n",
"from IPython.core.display import display, HTML\n",
"init_notebook_mode(connected=True)\n",
"import scipy.io as sio\n",
"import random\n",
"random.seed(123)\n",
"from plotly.subplots import make_subplots\n",
"\n",
"from repo2data.repo2data import Repo2Data\n",
"import os \n",
"data_req_path = os.path.join(\"..\", \"binder\", \"data_requirement.json\")\n",
"repo2data = Repo2Data(data_req_path)\n",
"data_path = repo2data.install()[0]\n",
"derivativesDir = os.path.join(data_path,\"qMRFlow\")\n",
"\n",
"sessions = ['vendor750rev','vendorPRIrev','vendorSKYrev',\n",
" 'rth750rev','rthPRIrev','rthSKYrev']\n",
"sides = [\"negative\",\"negative\",\"negative\",\"positive\",\"positive\",\"positive\"]\n",
"pointpos = [0,0,-0.14,0.14,0,0]\n",
"convert_tag = {'rth750':'G1VENUS','rthPRI':'S1VENUS','rthSKY':'S2VENUS',\n",
" 'vendor750':'G1NATIVE','vendorPRI':'S1NATIVE','vendorSKY':'S2NATIVE'}\n",
"colors = {'rth750':'rgba(255,8,30,1)','rthPRI':'rgba(255,117,0,1)','rthSKY':'rgba(255,220,0,1)',\n",
" 'vendor750':'rgba(0,75,255,1)','vendorPRI':'rgba(24,231,234,1)','vendorSKY':'rgba(24,234,141,1)'}\n",
"colors_line = {'rth750':'rgba(255,8,30,0.35)','rthPRI':'rgba(255,117,0,0.35)','rthSKY':'rgba(255,220,0,0.35)',\n",
" 'vendor750':'rgba(0,75,255,0.35)','vendorPRI':'rgba(24,231,234,0.35)','vendorSKY':'rgba(24,234,141,0.35)'}\n",
"\n",
"def plot_ridgelines(subject,sessions,metric,region,colors,sides,sampN,leg,lowlim,uplim):\n",
" df = pd.DataFrame()\n",
" traces = []\n",
" for jj in range(0,len(sessions),1):\n",
" avg_array = []\n",
" std_array = []\n",
" cur_region = derivativesDir + '/sub-' + subject + '/ses-' + sessions[jj] + '/stat/sub-' + subject + '_ses-' + sessions[jj] + '_desc-' + region + '_metrics.mat'\n",
" cur_contents = sio.loadmat(cur_region)\n",
" cur_contents = cur_contents[metric]\n",
" cur_contents = np.delete(cur_contents, np.where(cur_contents < lowlim))\n",
" cur_contents = np.delete(cur_contents, np.where(cur_contents > uplim))\n",
" cur_tag = sessions[jj][0:-3]\n",
" cur_contents = random.sample(list(cur_contents),sampN)\n",
" traces.append(go.Violin(\n",
" x=cur_contents,\n",
" name = convert_tag[cur_tag],\n",
" marker=dict(\n",
" opacity = 0.7, \n",
" outliercolor='rgba(255, 255, 255, 0.6)',\n",
" line=dict(\n",
" outliercolor='rgba(219, 64, 82, 0.6)',\n",
" outlierwidth=2)),\n",
" line_color= colors[cur_tag],\n",
" legendgroup= convert_tag[cur_tag],\n",
" showlegend=leg,\n",
" side=sides[jj],\n",
" pointpos = pointpos[jj]))\n",
" return traces\n",
"\n",
"def plot_mosaic(metric,axis_title,lowlim,uplim):\n",
" nSamp = 10000\n",
" traces1 = plot_ridgelines('invivo1',sessions,metric,'wm',colors,sides,nSamp,True,lowlim,uplim)\n",
" traces2 = plot_ridgelines('invivo1',sessions,metric,'gm',colors,sides,nSamp,False,lowlim,uplim)\n",
" traces3 = plot_ridgelines('invivo2',sessions,metric,'wm',colors,sides,nSamp,False,lowlim,uplim)\n",
" traces4 = plot_ridgelines('invivo2',sessions,metric,'gm',colors,sides,nSamp,False,lowlim,uplim)\n",
" traces5 = plot_ridgelines('invivo3',sessions,metric,'wm',colors,sides,nSamp,False,lowlim,uplim)\n",
" traces6 = plot_ridgelines('invivo3',sessions,metric,'gm',colors,sides,nSamp,False,lowlim,uplim)\n",
" fig = make_subplots(rows=3, cols=2,\n",
" subplot_titles=(\"P1 White Matter\", \"P1 Gray Matter\", \"P2 White Matter\", \"P2 Gray Matter\", \"P3 White Matter\", \"P3 Gray Matter\"),\n",
" vertical_spacing = 0.08,\n",
" shared_xaxes=True)\n",
" for trace in traces1:\n",
" fig.add_trace(trace,row=1, col=1)\n",
" for trace in traces2:\n",
" fig.add_trace(trace,row=1, col=2)\n",
" for trace in traces3:\n",
" fig.add_trace(trace,row=2, col=1)\n",
" for trace in traces4:\n",
" fig.add_trace(trace,row=2, col=2)\n",
" for trace in traces5:\n",
" fig.add_trace(trace,row=3, col=1)\n",
" for trace in traces6:\n",
" fig.add_trace(trace,row=3, col=2)\n",
" axis_template = dict(linecolor = 'white', \n",
" showticklabels = True,\n",
" tickfont=dict(color=\"white\",size=12), \n",
" gridcolor = 'rgba(255,255,255,0.2)')\n",
" fig.update_yaxes(axis_template,showgrid=False)\n",
" fig.update_traces(orientation='h')\n",
" fig.update_layout(paper_bgcolor = '#1c2a3d',plot_bgcolor = '#1c2a3d')\n",
" fig.update_traces(meanline=dict(visible=True,color=\"white\",width=3),\n",
" points= 'outliers', # show all points\n",
" jitter=0.1,\n",
" width=3) # add some jitter on points for better visibility\n",
" #scalemode='count')\n",
" fig.update_xaxes(axis_template)\n",
" fig.update_layout(width=800,height=950,title=f\"VENUS vs NATIVE in-vivo {metric} distributions for each participant (P)\", title_font_color=\"white\")\n",
" fig.update_annotations(font=dict(color=\"white\", family = \"Courier\"))\n",
" for i in range(1,7): \n",
" fig['layout']['xaxis{}'.format(i)]['title']= axis_title\n",
" fig.for_each_xaxis(lambda axis: axis.title.update(font=dict(color = 'white', size=12,family=\"Courier\"),standoff=0))\n",
" fig.update_layout(\n",
" legend=dict(\n",
" y = 0.5,\n",
" font=dict(\n",
" family=\"Courier\",\n",
" size=16,\n",
" color=\"white\")))\n",
" return fig\n",
"\n",
"import warnings\n",
"warnings.filterwarnings('ignore')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### VENUS vs NATIVE `T1` distributions\n",
"\n",
"`````{admonition} Note\n",
":class: tip\n",
"Reproduces Figure 5d-g from the [article](https://www.biorxiv.org/content/10.1101/2021.12.27.474259v2) for the Participant-3.\n",
"`````\n",
"\n",
"```{figure} ../assets/invivo_t1_p3.jpg\n",
"---\n",
"height: 400px\n",
"name: p3t1\n",
"---\n",
"vendor-native and VENUS quantitative T1 maps from P3 are shown in one axial slice.\n",
"```\n",
"\n",
"> The following KDEs agree well with the qualitative observations from the **Fig. 1** above: Inter-vendor agreement (G1-vs-S1 and G1-vs-S2) of VENUS is superior to that of vendor-native T1 maps, both in the GM and WM."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"fig = plot_mosaic('T1','T1 (s)',0.5,3)\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### VENUS vs NATIVE `MTR` distributions\n",
"\n",
"\n",
"`````{admonition} Note\n",
":class: tip\n",
"Reproduces Figure 5e-h from the [article](https://www.biorxiv.org/content/10.1101/2021.12.27.474259v2) for the Participant-3.\n",
"`````\n",
"\n",
"```{figure} ../assets/invivo_mtr_p3.jpg\n",
"---\n",
"height: 400px\n",
"name: p3mtr\n",
"---\n",
"nendor-native and VENUS quantitative MTR maps from P3 are shown in one axial slice.\n",
"```\n",
"\n",
"> The following KDEs agree well with the qualitative observations from the **Fig. 2** above: Inter-vendor agreement (G1-vs-S1 and G1-vs-S2) of VENUS is superior to that of vendor-native MTR maps, both in the GM and WM."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"fig = plot_mosaic('MTR','MTR (%)',35,70)\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### VENUS vs NATIVE `MTsat` distributions\n",
"\n",
"\n",
"`````{admonition} Note\n",
":class: tip\n",
"Reproduces Figure 5f-i from the [article](https://www.biorxiv.org/content/10.1101/2021.12.27.474259v2) for the Participant-3.\n",
"`````\n",
"\n",
"```{figure} ../assets/invivo_mtr_p3.jpg\n",
"---\n",
"height: 400px\n",
"name: p3mtsat\n",
"---\n",
"vendor-native and VENUS quantitative MTsat maps from P3 are shown in one axial slice.\n",
"```\n",
"\n",
"> The following KDEs agree well with the qualitative observations from the **Fig. 3** above: Inter-vendor agreement (G1-vs-S1 and G1-vs-S2) of VENUS is superior to that of vendor-native MTsat maps, both in the GM and WM."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"fig = plot_mosaic('MTsat','MTsat (a.u.)',0.5,5)\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Deciles and shift function\n",
"\n",
"\n",
"`````{admonition} Note\n",
":class: tip\n",
"Reproduces Figure 6a from the [article](https://www.biorxiv.org/content/10.1101/2021.12.27.474259v2).\n",
"`````\n",
"\n",
"Before moving on to the next page that performs shift function (for P3) and hierarchical shift function analyses (across participants) in `R`, this subsection provides some background information about these statistical tests.\n",
"\n",
"```{figure} ../assets/sf_exp.jpg\n",
"---\n",
"height: 600px\n",
"name: sfexp\n",
"---\n",
"explanation of the shift function analysis.\n",
"```\n",
"\n",
"### What is a decile? \n",
"\n",
"Deciles (a.k.a quantiles) are 1/10th of a distribution, created by splitting the distribution at 9 boundary points (excluding the min and max). For example, median is the 5th decile of a distribution.\n",
"\n",
"\n",
"### How do deciles create a shift function? \n",
"\n",
"To characterize differences between two distributions not only at one decile (e.g., the median), but throughout the whole range of the distributions, shift functions calculate differences at each decile. \n",
"\n",
"For example, distributions `X` (green) and `Y` (red) below do not differ much at the central tendency, marked by the dashed lines. The distribution `Y` is slightly higher than `X` at the central decile. \n",
"\n",
"On the other hand, distributions `X` and `Y` differ in spread. The distribution `X` is much wider than the distribution `Y`. As a result:\n",
"\n",
"* The first decile of `X` (green vertical line on the left) is smaller than that of `Y` (red vertical line on the left)\n",
"* The last decile of `X` (green vertical line on the right) is greater than that of `Y` (red vertical line on the right)\n",
"\n",
"The profile of shift function calculated for **X - Y** explains how the deciles of distribution `X` should be shifted to match those of the distribution `Y`: \n",
"\n",
"* The first decile of `Y` (**D1**) must be shifted towards left (towards lower values) to match that of `X`\n",
"* The first decile of `Y` must be shifted towards right (towards higher values) to match that of `X`\n",
"\n",
"Therefore, the first decile on Fig. 4 is below and the last decile is above the zero-crossing.\n",
"\n",
"`````{admonition} See Also\n",
":class: seealso\n",
"For a more detailed explanation on shift functions, please see this [blog post](https://garstats.wordpress.com/2016/07/12/shift-function/) by Guillaume A. Rousselet.\n",
"`````"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"scrolled": true,
"tags": [
"hide-input"
]
},
"outputs": [],
"source": [
"from scipy.stats import norm\n",
"# generate random numbers from N(0,1)\n",
"\n",
"fig = go.Figure()\n",
"\n",
"\n",
"data_normal = norm.rvs(size=10000,loc=0.5,scale=1)\n",
"data_normal2 = norm.rvs(size=10000,loc=0,scale=2)\n",
"\n",
"# Figure for ideal case \n",
"#data_normal = norm.rvs(size=10000,loc=0,scale=1)\n",
"#data_normal2 = norm.rvs(size=10000,loc=0,scale=1)\n",
"\n",
"fig.add_trace(go.Violin(x=data_normal,\n",
" opacity = 1, \n",
" line_color= \"crimson\",\n",
" side=\"positive\",\n",
" name = \"Y\",\n",
" meanline=dict(visible=True,color=\"crimson\",width=4),\n",
" width=2\n",
"))\n",
"fig.add_trace(go.Violin(x=data_normal2,\n",
" marker=dict(\n",
" opacity = 0.3, \n",
" outliercolor='rgba(255, 255, 255, 0.6)',\n",
" line=dict(\n",
" outliercolor='rgba(219, 64, 82, 0.6)',\n",
" outlierwidth=2)),\n",
" line_color= \"limegreen\",\n",
" side=\"positive\",\n",
" name = \"X\",\n",
" meanline=dict(visible=True,color=\"limegreen\",width=4),\n",
" width=2 \n",
"))\n",
"\n",
"\n",
"fig.update_traces(orientation='h')\n",
"fig.update_layout(paper_bgcolor = 'white',plot_bgcolor = 'white',height=600, width=800)\n",
"fig.update_traces(points=False)\n",
"qs = pd.qcut(data_normal2, 10, labels=False)\n",
"aa = []\n",
"for ii in range(10):\n",
" aa.append(np.median(data_normal2[qs==ii]))\n",
"fig.add_vrect(\n",
" x0=aa[9] + 0.8, x1=aa[9]+0.9,\n",
" fillcolor=\"limegreen\", opacity=0.7,\n",
" layer=\"above\", line_width=0,\n",
")\n",
"\n",
"for ii in range(0,9,1):\n",
" fig.add_vrect(\n",
" x0=aa[ii], x1=aa[ii] + 0.02,\n",
" fillcolor=\"limegreen\", opacity=0.3,\n",
" layer=\"above\", line_width=0)\n",
"fig.add_vrect(\n",
" x0=aa[0], x1=aa[0] + 0.1,\n",
" fillcolor=\"limegreen\", opacity=0.7,\n",
" layer=\"above\", line_width=0,\n",
")\n",
"qs = pd.qcut(data_normal, 10, labels=False)\n",
"aa = []\n",
"for ii in range(10):\n",
" aa.append(np.median(data_normal[qs==ii])) \n",
"fig.add_vrect(\n",
" x0=aa[9], x1=aa[9]+0.1,\n",
" fillcolor=\"crimson\", opacity=0.7,\n",
" layer=\"above\", line_width=0,\n",
")\n",
"fig.add_vrect(\n",
" x0=aa[0], x1=aa[0]+0.1,\n",
" fillcolor=\"crimson\", opacity=0.7,\n",
" layer=\"above\", line_width=0,\n",
")\n",
"for ii in range(0,9,1):\n",
" fig.add_vrect(\n",
" x0=aa[ii], x1=aa[ii] + 0.02,\n",
" fillcolor=\"crimson\", opacity=0.3,\n",
" layer=\"above\", line_width=0)\n",
" \n",
"fig.add_annotation(x=aa[0], y=0.8,\n",
" text=\"<-- D1\",\n",
" showarrow=False,\n",
" yshift=10)\n",
"fig.add_annotation(x=aa[-1], y=0.8,\n",
" text=\"D9-->\",\n",
" showarrow=False,\n",
" yshift=10)\n",
"\n",
"\n",
"\n",
"fig.show()"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### What about hierarchical shift functions? \n",
"\n",
"Hierarchical shift function (HSF) is an extension of the shift function analysis for the comparison of dependent distributions across multiple participants.\n",
"\n",
"In our case, each between-scanner comparison comes from a repeated measurement, given that the same participants were scanned using different scanners and pulse sequence implementations. \n",
"\n",
"`````{admonition} See Also\n",
":class: seealso\n",
"For a more detailed explanation on shift functions, please see this [blog post](https://garstats.wordpress.com/2019/02/21/hsf/) by Guillaume A. Rousselet.\n",
"`````"
]
}
],
"metadata": {
"celltoolbar": "Edit Metadata",
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.6.12"
}
},
"nbformat": 4,
"nbformat_minor": 4
}