# --- # jupyter: # jupytext: # formats: ipynb,py:percent # text_representation: # extension: .py # format_name: percent # format_version: '1.3' # jupytext_version: 1.13.6 # kernelspec: # display_name: Python 3 (ipykernel) # language: python # name: python3 # --- # %% [markdown] # ## Import Packages # %% # %matplotlib inline from typing import Callable, Iterable import nilearn from nilearn import plotting import nilearn.datasets as datasets import nilearn.image as image from nilearn import surface from sklearn.utils import shuffle import nibabel as nib import numpy as np import pandas as pd import seaborn as sns import matplotlib.pyplot as plt from matplotlib.colors import LinearSegmentedColormap from matplotlib.colors import ListedColormap from brainspace.gradient import GradientMaps from brainspace.utils.parcellation import map_to_labels from neurolang.frontend import NeurolangPDL, ExplicitVBR import data_utils from data_utils import vox2mm, merge_difumo_network_label, truncate_colormap, create_coordinate_dataset # %% [markdown] # ## Call NeuroLang's Probablistic Frontend # %% nl = NeurolangPDL() # %% [markdown] # ## Load Grey Matter Mask and Add it as Tuple List (i.e. Table) in NeuroLang # %% # Grey Matter Mask mni_mask = image.resample_img( nib.load(datasets.fetch_icbm152_2009()["gm"]), np.eye(3)*3, interpolation="nearest" ) plotting.plot_stat_map(mni_mask, title='Grey Matter Mask') # Extract GM voxel coordinates (probability > 0.25) and store them in table "GreyMatterVoxel" gm_matrix_coordinates = np.array(np.where(mni_mask.get_fdata() > 0.25)).T gm_xyz_coordinates = pd.DataFrame(vox2mm(gm_matrix_coordinates, affine=mni_mask.affine), columns=["x", "y", "z"]) nl.add_tuple_set(gm_xyz_coordinates, name="GreyMatterVoxel") # %% [markdown] # ## Load LPFC mask and Add it as a Table in NeuroLang # %% # Lateral Prefrontal Cortex Mask lpfc_mask=image.resample_to_img(source_img=nib.load('lpfc_mask.nii'), target_img=mni_mask, interpolation='nearest') plotting.plot_roi(lpfc_mask, cut_coords=(45,30,30), title='Lateral Prefrontal Cortex Mask') # Extract LPFC voxel coordinates and store them in table "LateralPFCVoxel" lpfc_matrix_coordinates=np.array(np.where(lpfc_mask.get_fdata() != 0)).T lpfc_xyz_coordinates = pd.DataFrame(vox2mm(lpfc_matrix_coordinates, affine=mni_mask.affine), columns = ["x", "y", "z"]) nl.add_tuple_set(lpfc_xyz_coordinates, name="LPFCVoxel") # %% [markdown] # ## Load Neurosynth Data and Add them to NeuroLang # %% # Neurosynth Data term_in_study, peak_reported, study_ids = data_utils.fetch_neurosynth( tfidf_threshold=0.001, ) # Store Peak Activation Coordinates in a "PeakReported" table nl.add_tuple_set(peak_reported, name="PeakReported") # Store Study IDs in a "Study" table nl.add_tuple_set(study_ids.drop_duplicates(), name="Study") # Add uniform weights for study selection in the meta-analysis and store them in "SelectedStudy" table nl.add_uniform_probabilistic_choice_over_set(study_ids.drop_duplicates(), name="SelectedStudy") # %% [markdown] # ## Load Difumo-1024 Atlas Data and Add them to NeuroLang # %% #Difumo 1024 Regions Atlas region_voxels, difumo_meta = data_utils.fetch_difumo( mask=mni_mask, coord_type="xyz", n_components=1024, ) nl.add_tuple_set(region_voxels.drop_duplicates(), name="RegionVoxel") # %% [markdown] # ## Define Functions to be used in NeuroLang Programs # %% EPS = 1e-20 @nl.add_symbol def log_odds(p, p0): a = min(p/(1-p) + EPS, 1.) b = min(p0/(1-p0) + EPS, 1.) logodds = np.log10(a/b) return logodds @nl.add_symbol def weighted_center(coords: Iterable[int], w: Iterable[float]) -> int: return int(np.sum(np.multiply(w, coords)) / np.sum(w)) @nl.add_symbol def agg_int_to_str(a,b): d = min(a, b) e = max(a,b) return str("%d-%d"%(d, e)) @nl.add_symbol def Euclidean_Dist(a, b, c, d, e, f): x = np.array((a, b, c)) y = np.array((d, e, f)) temp = x-y dist = np.sqrt(np.dot(temp.T, temp)) return dist @nl.add_symbol def maximum(w: Iterable[float]) -> float: return np.max(w) # %% [markdown] # ## Infer the VoxelReported table # %% with nl.scope as e: voxel_reported=nl.execute_datalog_program(""" VoxelReported(x, y, z, study) :- GreyMatterVoxel(x, y, z) & PeakReported(x2, y2, z2, study) & (distance == EUCLIDEAN(x, y, z, x2, y2, z2)) & (distance < 10) ans(x, y, z, study) :- VoxelReported(x, y, z, study) """) nl.add_tuple_set(voxel_reported.as_pandas_dataframe().drop_duplicates(), name="VoxelReported") # %% [markdown] # ## Finding Regions within the LPFC mask # %% with nl.scope as e: lpfc_regions = nl.execute_datalog_program(""" RegionVolume(r, count(x, y, z)) :- RegionVoxel(r, x, y, z, w) VolumeOfOverlapWithMask(r, count(x, y, z)) :- RegionVoxel(r, x, y, z, w) & LPFCVoxel(x, y, z) LPFCRegion(r) :- RegionVolume(r, v0) & VolumeOfOverlapWithMask(r, v) & (v/v0 > 0.5) ans(r) :- LPFCRegion(r) """) nl.add_tuple_set(lpfc_regions.as_pandas_dataframe().drop_duplicates(), name="LPFCRegion") # %% [markdown] # ## Inferring MetaAnalytic Connectivity Matrix # %% with nl.scope as e: connectivity_matrix = nl.execute_datalog_program(""" RegionMaxWeight(r, max(w)) :- RegionVoxel(r, x, y, z, w) RegionVoxelNormalizedWeight(r, x, y, z) :: w/W :- RegionVoxel(r, x, y, z, w) & RegionMaxWeight(r, W) LPFCRegionActive(r, study) :- RegionVoxelNormalizedWeight(r, x, y, z) & VoxelReported(x, y, z, study) & LPFCRegion(r) LPFCRegionNotActive(r, study) :- LPFCRegion(r) & ~LPFCRegionActive(r, study) & Study(study) BrainRegionActive(r, study) :- VoxelReported(x, y, z, study) & RegionVoxelNormalizedWeight(r, x, y, z) ProbabilityOfCoactivation(r, r2, PROB) :- BrainRegionActive(r, study) // (LPFCRegionActive(r2, study) & SelectedStudy(study)) ProbabilityOfNoCoactivation(r, r2, PROB) :- BrainRegionActive(r, study) // (LPFCRegionNotActive(r2, study) & SelectedStudy(study)) MetaAnalyticConnectivityMatrix(r2, r, max(p1, p0)) :- ProbabilityOfCoactivation(r, r2, p1) & ProbabilityOfNoCoactivation(r, r2, p0) ans(r2, r, LOR) :- MetaAnalyticConnectivityMatrix(r2, r, LOR) """) # %% [markdown] # ## Estimating the Gradients Using Diffusion Embedding # %% #compute similarity matrix using eta-squared metric and then derive gradients using PCA or DifussionEmbedding def eta2(X): S = np.zeros((X.shape[0],X.shape[0])) for i in range(0,X.shape[0]): for j in range(i,X.shape[0]): mi = np.mean([X[i,:],X[j,:]],0) mm = np.mean(mi) ssw = np.sum(np.square(X[i,:]-mi) + np.square(X[j,:]-mi)) sst = np.sum(np.square(X[i,:]-mm) + np.square(X[j,:]-mm)) S[i,j] = 1-ssw/sst S += S.T S -= np.eye(S.shape[0]) return S similarity_mat = eta2(connectivity_matrix.fillna(0).values) gradient_maps = GradientMaps(n_components=20, approach="dm") gradient_maps.fit(similarity_mat, sparsity=0.0) grads = pd.DataFrame(gm.gradients_, index=coactivation_mat.index) grads_join = grads.join(region_voxels.set_index('region')) i, j, k = np.round( nib.affines.apply_affine( np.linalg.inv(mni_mask.affine), grads_join[["x", "y", "z"]].values )).T.astype(int) gradient = np.zeros(mni_mask.shape) gradient[i, j, k] = pd.qcut(grads_join[0], 5, labels=np.arange(1, 6)) lpfc_mask = np.zeros_like(gradient) lpfc_mask[i, j, k] = 1 gradient = nib.Nifti1Image(gradient, affine=mni_mask.affine) lpfc_mask = nib.Nifti1Image(lpfc_mask, affine=mni_mask.affine) plotting.plot_img_on_surf(gradient, mask_img=lpfc_mask, threshold=1, symmetric=False, cmap="gnuplot") plt.show() plt.plot(gm.lambdas_ / gm.lambdas_.sum()); plt.show()