TLE-HC

Overview

Case-control difference mapping in temporal lobe epilepsy. This notebook will examine morphological features of the hippocampus in TLE patients. We compare the ipsilateral foci hippocampi to contralateral (asymmetry analysis), and to age-matched controls (case-control analysis).

import numpy as np
import matplotlib.pyplot as plt
import nibabel as nib
import hippomaps as hm
import pandas as pd
from scipy.stats import ttest_1samp
import time
start_time = time.time()

# config

useCheckpoints = True
if useCheckpoints:
    hm.fetcher.get_tutorialCheckpoints(['TLE-HC.npz', 'participants_MICs_PX.csv', 'participants_MICs_HC.csv'])

hippunfold_dir = '/data/mica3/BIDS_MICs/derivatives/hippunfold_v1.3.0/hippunfold'

subs_PX = pd.read_csv(f'checkpoints/participants_MICs_PX.csv', delimiter=',')
subs_HC = pd.read_csv(f'checkpoints/participants_MICs_HC.csv', delimiter=',')
hemis = ['L','R']
labels = ['hipp']#,'dentate']
den = '0p5mm'

# get expected number of vertices and their indices
nV,iV = hm.config.get_nVertices(labels,den)

# loop through these types of structural images
features = ['thickness','gyrification','curvature','surfarea']

0) Map and load volumetric data to surfaces

As in all tutorials here, this step is OPTIONAL, and will be skipped by default. It provides an example of how data can be mapped to hippocampal surfaces outside of python (using wb_command). This relies on having the data stored locally, and should be considered example code. This code may differ depending on where/how your data is stored and formatted, and so may require some customization for new projects. For the purposes of this tutorial, we provide a matrix of loaded data at the end, so skip to the next step.

In this example, we loop through samples (that is, subjects and hemipsheres) pulling hippunfold output features (that is, features intrinsic to the surfaces themselves, such as thickness).

if not useCheckpoints:

    # Load PX data
    hipp_dat_PX = np.zeros([nV,len(hemis),len(subs_PX['ID']),len(features)])*np.nan
    for f,feature in enumerate(features):
        for s,sub in enumerate(subs_PX['ID']):
            for h,hemi in enumerate(hemis):
                for l,label in enumerate(labels):
                    ses = subs_PX['SES'][s]
                    if subs_PX['pathology'][s][0] == hemi:
                        ic = 0 # always put ipsi first
                    else:
                        ic=1 # contra second
                    fn = f'{hippunfold_dir}/{sub}/{ses}/surf/{sub}_{ses}_hemi-{hemi}_space-T1w_den-{den}_label-{label}_{feature}.shape.gii'
                    if feature=='curvature' and hemi == 'R':
                        hipp_dat_PX[iV[l],ic,s,f] = -nib.load(fn).darrays[0].data
                    else:
                        hipp_dat_PX[iV[l],ic,s,f] = nib.load(fn).darrays[0].data

    # Load HC data
    hipp_dat_HC = np.zeros([nV,len(hemis),len(subs_HC['ID']),len(features)])*np.nan
    for f,feature in enumerate(features):
        for s,sub in enumerate(subs_HC['ID']):
            for h,hemi in enumerate(hemis):
                for l,label in enumerate(labels):
                    ses = subs_HC['SES'][s]
                    fn = f'{hippunfold_dir}/{sub}/{ses}/surf/{sub}_{ses}_hemi-{hemi}_space-T1w_den-{den}_label-{label}_{feature}.shape.gii'
                    if feature=='curvature' and hemi == 'R':
                        hipp_dat_HC[iV[l],h,s,f] = -nib.load(fn).darrays[0].data
                    else:
                        hipp_dat_HC[iV[l],h,s,f] = nib.load(fn).darrays[0].data

    np.savez_compressed("checkpoints/TLE-HC",hipp_dat_PX, hipp_dat_HC)

# overview demographics of patients

print(f"mean age: {np.mean(subs_PX['AGE'])} +/- {np.std(subs_PX['AGE'])}")
print(f"M/F: {np.sum(subs_PX['SEX']=='M')}/{np.sum(subs_PX['SEX']=='F')}")
print(f"L/R: {np.sum(subs_PX['pathology']=='Right TLE')}/{np.sum(subs_PX['pathology']=='Left TLE')}")

mean age: 42.43478260869565 +/- 11.021031621233071
M/F: 12/11
L/R: 10/13

# overview demographics of controls

print(f"mean age: {np.mean(subs_HC['AGE'])} +/- {np.std(subs_HC['AGE'])}")
print(f"M/F: {np.sum(subs_HC['SEX']=='M')}/{np.sum(subs_HC['SEX']=='F')}")

mean age: 40.42857142857143 +/- 10.43378864929855
M/F: 21/21

# Here we calculate and add a new feature: columnar volume. This is simply thicknes (mm) * surface area (mm2) giving a units in mm3.

hipp_dat_PX = np.load("checkpoints/TLE-HC.npz")['arr_0']
hipp_dat_HC= np.load("checkpoints/TLE-HC.npz")['arr_1']

coolumnar_volume = hipp_dat_PX[:,:,:,0]*hipp_dat_PX[:,:,:,3]
hipp_dat_PX = np.concatenate((hipp_dat_PX,coolumnar_volume.reshape([nV,len(hemis),len(subs_PX['ID']),1])),axis=3)

coolumnar_volume = hipp_dat_HC[:,:,:,0]*hipp_dat_HC[:,:,:,3]
hipp_dat_HC = np.concatenate((hipp_dat_HC,coolumnar_volume.reshape([nV,len(hemis),len(subs_HC['ID']),1])),axis=3)

features.append('columnarVol')

1) plot group-averages

Here we plot averages of the ipsilateral seizure focus, contralalteral, and healthy controls. We then compare ipsilateral to contralalteral (an asymmetry analysis) and ipsilateral to controls (case-control differences). Asymmetry is nice since it includes perfect matching on things like demographics, etc., but in cases with severe epilepsy there can also be abnormalities or secondary epileptogenic zones in the contralateral hemisphere as well. Case-control differences are less closely matched, but we can be confident that the HC data is normative.

ipsilateral = np.mean(hipp_dat_PX[:,0,:,:], axis=1)
contralateral = np.mean(hipp_dat_PX[:,1,:,:], axis=1)
ctrl = np.mean(hipp_dat_HC,axis=(1,2))
asymmetry = contralateral-ipsilateral
caseCtrl = ctrl-ipsilateral

cdata1 = np.stack((ipsilateral,contralateral,ctrl),axis=2) # these will share color range
cdata2 = np.stack((asymmetry,caseCtrl),axis=2) # these will share their own separate color range

# plot altogether
fig, ax = plt.subplots(2,5, figsize=(25,25))
for f in range(5):
    hm.plotting.surfplot_canonical_foldunfold(cdata1[:,f,:], color_bar=('right'), hemis=['L'], labels=['hipp'], cmap='Blues', color_range='sym', unfoldAPrescale=True, share='both', tighten_cwindow=True, embed_nb=True, screenshot=True, filename='tmp.png')
    i = plt.imread('tmp.png')
    ax[0,f].imshow(i)
    ax[0,f].set_axis_off()
    ax[0,f].set_anchor("NW")
    ax[0,f].set_title(features[f])
    hm.plotting.surfplot_canonical_foldunfold(cdata2[:,f,:], color_bar=('right'), hemis=['L'], labels=['hipp'], cmap='bwr', color_range='sym', unfoldAPrescale=True, share='both', tighten_cwindow=True, embed_nb=True, screenshot=True, filename='tmp.png')
    i = plt.imread('tmp.png')
    ax[1,f].imshow(i)
    ax[1,f].set_axis_off()
    ax[1,f].set_anchor("NW")

!rm tmp.png

2) Examine inter-patient (ipsilateral) consistency

Here, we focus only on case-control differences for each ipsilateral hippocampus to the control average.

# make a histogram of the off-diagonal cross-patient correlations (minus controls average)
mfcorr = []
sdfcorr = []
allcorr = []
fig, ax = plt.subplots(nrows=1, ncols=len(features), figsize=(3*len(features),3))
for f,feature in enumerate(features):
    cdat = (ctrl[:,f]-hipp_dat_PX[:,0,:,f].T).T
    corr = np.corrcoef(cdat.T) # as above 2
    fcorr = corr[np.triu_indices(20,k=1)] # this return only the off-diagonal (lower left triangle)
    print(ttest_1samp(fcorr,0,nan_policy='omit'))
    ax[f].hist(fcorr) # make a histogram for this particular feature
    mfcorr.append(np.mean(fcorr))
    sdfcorr.append(np.std(fcorr))
    allcorr.append(fcorr)

TtestResult(statistic=11.806947813278736, pvalue=1.86729012634132e-24, df=189)
TtestResult(statistic=15.077496927036245, pvalue=3.02929417496767e-34, df=189)
TtestResult(statistic=18.179538005866, pvalue=2.3078375623474734e-43, df=189)
TtestResult(statistic=17.15560123203877, pvalue=2.1642355212648486e-40, df=189)
TtestResult(statistic=16.024709512021285, pvalue=4.612713044556611e-37, df=189)

# Generate individual points from provided allcorr data
jitter_strength = 0.1  # Controls horizontal spread

data_points = []
x_positions = np.arange(len(features))
for i, points in enumerate(allcorr):
    jitter = np.random.uniform(-jitter_strength, jitter_strength, size=len(points))
    data_points.append((x_positions[i] + jitter, points))

# Plot
fig, ax = plt.subplots(figsize=(2 * len(features), 2))
for x, y in data_points:
    ax.scatter(x, y, color='black', alpha=0.2, s=10)  # Individual points in greyscale
ax.errorbar(x_positions, mfcorr, yerr=sdfcorr, fmt='o', color='black', capsize=5)  # Mean and SD line in greyscale

# Remove outer border
ax.spines['top'].set_visible(False)
ax.spines['right'].set_visible(False)
ax.spines['left'].set_visible(False)
ax.spines['bottom'].set_visible(False)

# Add horizontal line at 0
ax.axhline(0, color='black', linestyle='--', linewidth=1)

ax.set_xticks(x_positions)
ax.set_xticklabels(features)
plt.show()

3) Contextualize with other maps

t, ax = hm.stats.contextualize2D(caseCtrl[:,:2], taskNames=features[:2], numerbMaps=True, nperm=10000) # set nperm low for fast testing!

/export03/data/opt/venv/lib/python3.9/site-packages/numpy/lib/function_base.py:2897: RuntimeWarning: invalid value encountered in divide
  c /= stddev[:, None]
/export03/data/opt/venv/lib/python3.9/site-packages/numpy/lib/function_base.py:2898: RuntimeWarning: invalid value encountered in divide
  c /= stddev[None, :]

save

# save a copy of the 2D map
!mkdir -p ../maps/HippoMaps-initializationMaps/Dataset-MICs
for f,feature in enumerate(features):
    for h,hemi in enumerate(hemis):
        for l,label in enumerate(labels):
            cdat = caseCtrl[:,f]
            data_array = nib.gifti.GiftiDataArray(data=cdat.astype(np.float32))
            image = nib.gifti.GiftiImage()
            image.add_gifti_data_array(data_array)
            nib.save(image, f'../maps/HippoMaps-initializationMaps/Dataset-MICs/MRI-3T-MRI-TLE-case-control-{feature}_average-{len(subs_PX)}_hemi-mix_den-2mm_label-{label}.shape.gii')

end_time = time.time()
duration = end_time - start_time
print(f"Total duration: {duration:.2f} seconds")
duration

Total duration: 715.16 seconds

715.1646678447723