Neuro ML 2020

pip install nilearn

pip install diagram2vec

import os
import re

import pandas as pd
import numpy as np

import nibabel as nib

import nilearn
from nilearn.image import concat_imgs
from nilearn.input_data import NiftiLabelsMasker
from nilearn.image import high_variance_confounds
from nilearn.connectome import ConnectivityMeasure
from nilearn import datasets

%matplotlib inline
import matplotlib.pyplot as plt

import numpy as np
import networkx as nx 
import nilearn # pip install nilearn
from scipy.stats import spearmanr
from sklearn.metrics import mutual_info_score
import warnings
warnings.filterwarnings("ignore")

from google.colab import drive
drive.mount('/content/drive')

! unzip drive/MyDrive/нейро/neuro2022/seminar_5/abide2.zip && ls

data = 'ABIDEII-GU_1'

def make_correlation_matrix(path_to_fmriprep_data, path_to_save_connectivity_matrices=False):
    """
        Process the fmriprep preprocessed functional MRI time-series into 2D correlation matrix as DataFrame using Nilearn lib.
    """
   
    dataset = datasets.fetch_atlas_aal(version='SPM12', data_dir=None, url=None, resume=True, verbose=1)
    atlas_filename = dataset.maps
    labels = dataset.labels
    correlation_measure = ConnectivityMeasure(kind='correlation')
    ts =  []
    for n_patient, patient_folder_name in enumerate(os.listdir(path_to_fmriprep_data)[:4]):
        print(n_patient, patient_folder_name)
        temp_path = os.path.join(path_to_fmriprep_data, patient_folder_name)
        print(temp_path)
           # temp_path = path_to_fmriprep_data + patient_folder_name + '/func/'
        for filename in os.listdir(os.path.join(temp_path, 'ses-1/func')):
                    print(filename)

                    print(os.path.join(temp_path, 'ses-1/func', filename))

#                     tr = tr_extractor(os.path.join(temp_path, 'ses-1/func', filename))
#                     print(tr)
                    tr = 2
                     # extracted data from non-overlapping volumes
                    masker = NiftiLabelsMasker(labels_img=atlas_filename, background_label=labels, 
                               standardize=True, detrend=True, resampling_target='labels', 
                               low_pass=0.08, high_pass=0.009, t_r=tr, memory='nilearn_cache', memory_level=1, verbose=0)
                    
                    img = concat_imgs(os.path.join(temp_path, 'ses-1/func', filename), auto_resample=True, verbose=0)
                    print(img.shape)
                    confounds = high_variance_confounds(img, 1)
                    time_series = masker.fit_transform(img, confounds)
                    print(time_series.shape)
                    ts.append(time_series)
                    correlation_matrix = correlation_measure.fit_transform([time_series])[0]
                    np.fill_diagonal(correlation_matrix, 1)

    return ts, time_series, correlation_matrix    

# load ABIDE1 data from NYU site
ts, time_series, correlation_matrix = make_correlation_matrix(data)

ts = np.stack(ts)
(n_patients, n_steps, n_regions) = ts.shape

n_patients, n_steps, n_regions

# visualization of raw time series
plt.figure(figsize=(16.5,5))
plt.title("First 3 time series")
plt.hlines(0, -2, n_steps+2, linewidth=1.0, linestyles="dotted")
plt.plot(ts[0,:,:3])
plt.show()

# mean and standard deviation
print("Mean ± std of the:")
print("1st time series: {:.3f} ± {:.3f}".format(np.mean(ts[0,:,0]), np.std(ts[0,:,0])))
print("2nd time series: {:.3f} ± {:.3f}".format(np.mean(ts[0,:,1]), np.std(ts[0,:,1])))
print("3rd time series: {:.3f} ± {:.3f}".format(np.mean(ts[0,:,2]), np.std(ts[0,:,2])))

# break 2nd and 3rd time series
ts[0,:,1] = ts[0,:,1] + 150 # add mean shift
ts[0,:,2] = ts[0,:,2] + np.linspace(0, 150, n_steps) # add trend

print("Mean ± std of the:")
print("1st time series: {:.3f} ± {:.3f}".format(np.mean(ts[0,:,0]), np.std(ts[0,:,0])))
print("2nd time series: {:.3f} ± {:.3f}".format(np.mean(ts[0,:,1]), np.std(ts[0,:,1])))
print("3rd time series: {:.3f} ± {:.3f}".format(np.mean(ts[0,:,2]), np.std(ts[0,:,2])))

# visualization of broken time series
plt.figure(figsize=(16.5,5))
plt.title("First 3 time series")
plt.hlines(0, -2, n_steps+2, linewidth=1.0, linestyles="dotted")
plt.plot(ts[0,:,:3])
plt.show()

from nilearn.signal import clean

ts_normalized = np.zeros_like(ts)

# normalize and detrend
for i in range(ts.shape[0]):
    ts_normalized[i] = clean(ts[i], standardize="zscore", detrend=True)
    
print("Mean ± std of the:")
print("1st time series: {:.3f} ± {:.3f}".format(np.mean(ts_normalized[0,:,0]), np.std(ts_normalized[0,:,0])))
print("2nd time series: {:.3f} ± {:.3f}".format(np.mean(ts_normalized[0,:,1]), np.std(ts_normalized[0,:,1])))
print("3rd time series: {:.3f} ± {:.3f}".format(np.mean(ts_normalized[0,:,2]), np.std(ts_normalized[0,:,2])))

# visualization of raw time series
plt.figure(figsize=(16.5,5))
plt.title("First 3 time series")
plt.hlines(0, -2, n_steps+2, linewidth=1.0, linestyles="dotted")
plt.plot(ts_normalized[0,:,:3])
plt.show()

ts = ts_normalized

from nilearn.connectome import ConnectivityMeasure
from sklearn.covariance import EmpiricalCovariance

covariance_estimator = EmpiricalCovariance()
connectivity_correlation = ConnectivityMeasure(kind="correlation", cov_estimator=covariance_estimator)

ts[0].shape

R = connectivity_correlation.fit_transform(ts)
R.shape

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16.5,5))
plt.suptitle("Correlation connectivity matrices")
ax1.imshow(R[0])
ax2.imshow(R[1])
ax3.imshow(R[2])
plt.show()

np.min(np.linalg.eigvalsh(R[0]))

### your code here

alpha = 0.01
R_regularized = R[0] + alpha * np.eye(n_regions)

from sklearn.covariance import LedoitWolf

cov_estimator_shrinked = LedoitWolf()
connectivity_correlation_shrinked = ConnectivityMeasure(kind="correlation", cov_estimator=cov_estimator_shrinked)

# check the value of beta parameter
cov = cov_estimator_shrinked.fit(ts[0])
cov.shrinkage_

R_shrinked = connectivity_correlation_shrinked.fit_transform(ts)

# checking the minimum eigenvalue
np.min(np.linalg.eigvalsh(R_shrinked[0]))

from scipy.stats import spearmanr

S = np.zeros((3, n_regions, n_regions))

for i in range(3):
    S[i], _ = spearmanr(ts[i])

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16.5,5))
plt.suptitle("Spearman correlation connectivity matrices")
ax1.imshow(S[0])
ax2.imshow(S[1])
ax3.imshow(S[2])
plt.show()

plt.figure(figsize=(6,6))
plt.plot(ts[2,:,0], ts[2,:,12], ".")

from sklearn.metrics import mutual_info_score

def calc_MI(x, y, bins=10):
    c_xy = np.histogram2d(x, y, bins)[0]
    mi = mutual_info_score(None, None, contingency=c_xy)
    return mi

def bound(x):
    return np.sqrt(1 - np.exp(-2 * x))

l = calc_MI(ts[1,:,42], ts[1,:,45])
l, bound(l)

c_xy = np.rot90(np.histogram2d(ts[2,:,0], ts[2,:,12], 10)[0])
plt.imshow(c_xy)

l = calc_MI(ts[1,:,42], ts[1,:,41])
l, bound(l)

plt.figure(figsize=(6,6))
plt.plot(ts[2,:,42], ts[2,:,41], ".")

c_xy = np.rot90(np.histogram2d(ts[2,:,42], ts[2,:,41], 10)[0])
plt.imshow(c_xy)

%%time
M = np.zeros((3, n_regions, n_regions))

for k in range(3):
    for i in range(n_regions):
        for j in range(i, n_regions):
            M[k,i,j] = bound(calc_MI(ts[k,:,i], ts[k,:,j]))
            
    M[k] = M[k] + M[k].T
    np.fill_diagonal(M[k], 1)

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16.5,5))
plt.suptitle("Mutual information connectivity matrices")
ax1.imshow(M[0])
ax2.imshow(M[1])
ax3.imshow(M[2])
plt.show()

R_tresholded = R[0].copy()
np.fill_diagonal(R_tresholded, 0)
R_tresholded[R_tresholded < 0.5] = 0.0
R_tresholded

fig, (ax1) = plt.subplots(1, 1, figsize=(5,5))
plt.suptitle("Thresholded Pearson correlation matrix")
ax1.imshow(R_tresholded)
plt.show()

from nilearn import plotting

# get coordinates of brain regions
atlas_aal = nilearn.datasets.fetch_atlas_aal()
coordinates = plotting.find_parcellation_cut_coords(labels_img=atlas_aal["maps"])

coordinates.shape

# nilearn graph drawing
fig = plt.figure(figsize=(13,6))
edge_options = {"color": "r", "linewidth": 1.5, "alpha": 0.5}
plotting.plot_connectome(R_tresholded, coordinates, figure=fig, edge_kwargs=edge_options)

# matplotlib graph drawing
fig, (ax1) = plt.subplots(1, 1, figsize=(8,8))
ax1.set_title("Correlation graph")
nx.draw_shell(nx.from_numpy_array(R_tresholded), ax=ax1)
plt.show()

# create graph from connectivity matrix
G_R = nx.from_numpy_array(R_tresholded)

degree = np.array([degree[1] for degree in nx.degree(G_R)])
degree

neighbor_degree_avg = np.array(list(nx.average_neighbor_degree(G_R).values()))
neighbor_degree_avg

centrality_betweenness = np.array(list(nx.betweenness_centrality(G_R).values()))
centrality_betweenness

centrality_closeness = np.array(list(nx.closeness_centrality(G_R).values()))
centrality_closeness

clustering_coefficient_local = np.array(list(nx.clustering(G_R).values()))
clustering_coefficient_local

nx.local_efficiency(G_R)

nx.global_efficiency(G_R)

eigenvalues, _ = np.linalg.eigh(R[0])
eigenvalues

R_tresholded

# laplacian matrix L = D - A
# A is contained in variable R_tresholded
A = R_tresholded

# D, your code here
D = np.diag(A.sum(axis=0))

# L, your code here
L = D - A

eigenvalues_laplacian, _ = np.linalg.eigh(L)
eigenvalues_laplacian

plt.plot(eigenvalues_laplacian)
plt.show()

pip install ripser

from ripser import ripser
import diagram2vec

# reverse matrix to add higher correlated edges first to the filtration
R_filtered = 1 - np.abs(R[0])

# compute persistence diagram of the network
diagram_R = ripser(R_filtered, distance_matrix=True)["dgms"]
diagram_R

# vectorize persistent diagram
betti_curve = diagram2vec.persistence_curve(diagram_R)
betti_curve

fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(16.5,5))
ax1.set_title("Persistence diagram")
ax1.set_xlim(-0.025,1)
ax1.set_ylim(-0.025,1)
ax1.scatter(diagram_R[0][:,0], diagram_R[0][:,1], c="b")
ax1.scatter(diagram_R[1][:,0], diagram_R[1][:,1], c="r")
ax2.set_title("0th Betti number curve")
ax2.plot(betti_curve[0,0], c="b")
ax3.set_title("1st Betti number curve")
ax3.plot(betti_curve[0,1], c="r")
plt.show()

from sklearn.tree import DecisionTreeClassifier
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import cross_val_score
from tqdm import tqdm

def make_correlation_matrix(path_to_fmriprep_data, path_to_save_connectivity_matrices=False):
    """
        Process the fmriprep preprocessed functional MRI time-series into 2D correlation matrix as DataFrame using Nilearn lib.
    """
   
    dataset = datasets.fetch_atlas_aal(version='SPM12', data_dir=None, url=None, resume=True, verbose=1)
    atlas_filename = dataset.maps
    labels = dataset.labels
    correlation_measure = ConnectivityMeasure(kind='correlation')
    ts =  []
    for n_patient, patient_folder_name in tqdm(enumerate(os.listdir(path_to_fmriprep_data))):
        print(n_patient, patient_folder_name)
        temp_path = os.path.join(path_to_fmriprep_data, patient_folder_name)
        print(temp_path)
           # temp_path = path_to_fmriprep_data + patient_folder_name + '/func/'
        for filename in os.listdir(os.path.join(temp_path, 'ses-1/func')):
                    print(filename)

                    print(os.path.join(temp_path, 'ses-1/func', filename))

#                     tr = tr_extractor(os.path.join(temp_path, 'ses-1/func', filename))
#                     print(tr)
                    tr = 2
                     # extracted data from non-overlapping volumes
                    masker = NiftiLabelsMasker(labels_img=atlas_filename, background_label=labels, 
                               standardize=True, detrend=True, resampling_target='labels', 
                               low_pass=0.08, high_pass=0.009, t_r=tr, memory='nilearn_cache', memory_level=1, verbose=0)
                    
                    img = concat_imgs(os.path.join(temp_path, 'ses-1/func', filename), auto_resample=True, verbose=0)
                    print(img.shape)
                    confounds = high_variance_confounds(img, 1)
                    time_series = masker.fit_transform(img, confounds)
                    print(time_series.shape)
                    ts.append(time_series)
                    correlation_matrix = correlation_measure.fit_transform([time_series])[0]
                    np.fill_diagonal(correlation_matrix, 1)

    return ts, time_series, correlation_matrix    

# load ABIDE1 data from NYU site
ts, time_series, correlation_matrix = make_correlation_matrix(data)

ts = np.stack(ts)
(n_patients, n_steps, n_regions) = ts.shape

# compute correlation networks
R = ConnectivityMeasure(kind="correlation").fit_transform(ts)
R.shape

R.shape

# topological features
R_filtration = 1 - np.abs(R)

diagrams = []

for i, R_filtered in enumerate(R_filtration):
    diagram = ripser(R_filtered, distance_matrix=True)["dgms"]
    diagrams.append(diagram)

X_topological = diagram2vec.persistence_curve(diagrams, quantity="persistence", m=40)

X_topological[:,0].shape

dataset = datasets.fetch_atlas_aal(version='SPM12', data_dir=None, url=None, resume=True, verbose=1)
len(dataset.labels)

clf = DecisionTreeClassifier(random_state=42, max_depth=5)
# cross_val_score(clf, X_topological[:,1], dataset.labels[:16], cv=10).mean()

R_thresholded = np.copy(R)
for i in range(R_thresholded.shape[0]):
    np.fill_diagonal(R_thresholded[i], 0)

R_thresholded[R_thresholded < 0.47] = 0.0

# graph features
X_graph = np.zeros((n_patients, n_regions*5))

for i in range(R_thresholded.shape[0]):
    
    G = nx.from_numpy_array(R_thresholded[i])
    
    degree = np.array([degree[1] for degree in nx.degree(G)])
    neighbor_degree_avg = np.array(list(nx.average_neighbor_degree(G).values()))
    centrality_betweenness = np.array(list(nx.betweenness_centrality(G).values()))
    centrality_closeness = np.array(list(nx.closeness_centrality(G).values()))
    clustering_coefficient = np.array(list(nx.clustering(G).values()))
    
    feature_i = np.concatenate((degree, neighbor_degree_avg, centrality_betweenness, centrality_closeness, clustering_coefficient))
    X_graph[i] = feature_i

# cross_val_score(clf, X_graph, dataset.labels[:16], cv=10).mean()

# spectral features (eigher spectrum or Laplacian spectrum)
X_spectral = np.zeros((n_patients, n_regions))

for i in range(R_thresholded.shape[0]):
    eigenvalues_i, _ = np.linalg.eigh(R_thresholded[i])
    X_spectral[i] = eigenvalues_i

# cross_val_score(clf, X_spectral, dataset.labels[:16], cv=10).mean()

# cross_val_score(clf, np.concatenate((X_graph, X_spectral, X_topological[:,1]), axis=1), dataset.labels[:16], cv=10).mean()