Merge from master to tempFromCommit

ridge_cpm.py

@@ -0,0 +1,96 @@
+'''
+@Author: Or Duek
+Using Ridge regression to run the CPM for the KPE study
+'''
+#%% Load libraries
+from sklearn.linear_model import Ridge
+from sklearn.model_selection import GridSearchCV
+from sklearn.model_selection import cross_val_predict
+from sklearn.model_selection import cross_val_score
+from sklearn import linear_model
+from sklearn.model_selection import KFold
+from sklearn.model_selection import RepeatedKFold
+from sklearn.model_selection import PredefinedSplit
+from sklearn.feature_selection import SelectPercentile
+from sklearn.feature_selection import f_regression
+from sklearn.pipeline import Pipeline
+from sklearn.model_selection import StratifiedKFold
+from sklearn import metrics
+import pandas as pd
+import numpy as np
+import scipy.io as sio
+import h5py
+
+import time
+
+
+#%% Load matrices and behaviour
+# AAL
+first = np.load('/home/or/kpe_task_analysis/trauma_ses1.npy')
+second = np.load('/home/or/kpe_task_analysis/trauma_ses2.npy')
+delta = np.subtract(second, first)
+
+# delta1 = np.delete(delta, mask, axis=2)
+# delta1.shape
+y = np.array([-2,  1,  -30,  -17,  -28,   -4,  -30,  -18,  -22,  -18,  1,  -11,   -2,   -4, 16,  -32,   -8,  -14,  -20,   -3,  -23])
+delta.shape
+
+#%% turn matrix to array
+vecs = []
+for i in range(delta.shape[2]):
+    mat = delta[:,:,i]
+    flat = mat.flatten()
+    vecs.append(flat) 
+vecs = np.array(vecs)
+print(vecs.shape)
+vecs_reshape = np.moveaxis(vecs,0,-1)
+#%% Set parameters of the model
+pct = 0.8 # percent of edges kept in feature selection
+alphas = 10**np.linspace(10,-2,100)*0.5 # specify alphas to search
+#%%
+rg_grid = GridSearchCV(Ridge(normalize=False), cv=10, param_grid={'alpha':alphas}, iid=False)
+# using LASSO regression instead of ridge
+lasso = linear_model.Lasso
+lasso_grid = GridSearchCV(lasso(normalize=False), cv=10, param_grid={'alpha':alphas}, iid=False)
+
+reg = Pipeline([
+  ('feature_selection', SelectPercentile(f_regression, percentile=pct)),
+  ('regression', lasso_grid)
+])
+
+cv10 = KFold(n_splits=21)#, random_state=665)
+rpcv10 = RepeatedKFold(n_splits=3,n_repeats=3, random_state=665)
+# %% Run model
+start = time.time() # time the function
+all_pred = cross_val_predict(reg, vecs_reshape.T, y, cv=cv10, n_jobs=4)
+#all_score = cross_val_score(reg, vecs_reshape.T, y, cv=rpcv10, n_jobs=1) # repeated kfolds
+end = time.time()
+print(end - start) # print function running time
+
+# %%
+print(np.corrcoef(all_pred.T, y.T))
+
+# %%
+import scipy
+scipy.stats.pearsonr(all_pred, y)
+#%% plot
+import seaborn as sns
+import matplotlib.pyplot as plot
+sns.regplot(all_pred, y)
+
+
+#%% Try with shen parcellation
+first = np.load('/media/Data/Lab_Projects/KPE_PTSD_Project/neuroimaging/KPE_results/shen/trauma_ses1_shen.npy')
+second = np.load('/media/Data/Lab_Projects/KPE_PTSD_Project/neuroimaging/KPE_results/shen/trauma_ses2_shen.npy')
+delta = np.subtract(second, first)
+
+
+# %%
+# lets plot with groups also
+df = pd.DataFrame({'group': group_label, 'observed':y, 'predicted':all_pred})
+
+
+sns.lmplot('predicted','observed', hue='group', data= df)
+print(f'Ketamine group correlation {scipy.stats.pearsonr(df.predicted[df.group==1], df.observed[df.group==1])}')
+print(f'Midazolam group correlation {scipy.stats.pearsonr(df.predicted[df.group==0], df.observed[df.group==0])}')
+# %%

rs/.ipynb_checkpoints/difumo_extract_timeline-checkpoint.py

@@ -0,0 +1,112 @@
+# %%
+'''
+@author: Or Duek
+@date: Jul 16 2020
+
+This is a sciprt that uses the nee DiFuMo dictionary
+atlas (https://www.sciencedirect.com/science/article/pii/S1053811920306121#appsec7)
+
+In this file we will create a task based 
+'''
+# %% import libraries
+import pandas as pd 
+from nilearn.input_data import NiftiMapsMasker
+from nilearn import connectome
+from nilearn import datasets
+import numpy as np
+import nilearn.plotting
+from sklearn.model_selection import StratifiedShuffleSplit
+import os
+import glob
+from nilearn import connectome
+import seaborn as sns
+# %% Set output folder
+output_dir = '/gpfs/gibbs/pi/levy_ifat/Or/kpe/results/resting'
+# set session
+ses= '3' # session is a string
+# %% Functions
+# extract RS data and create vector for each subject
+def removeVars (confoundFile):
+    # this method takes the csv regressors file (from fmriPrep) and chooses a few to confound. You can change those few
+    import pandas as pd
+    confound = pd.read_csv(confoundFile,sep="\t", na_values="n/a")
+    finalConf = confound[['csf', 'white_matter', 'framewise_displacement', 
+                          'a_comp_cor_00', 'a_comp_cor_01','a_comp_cor_02', 'a_comp_cor_03', 
+                          'a_comp_cor_04', 'a_comp_cor_05', 'trans_x', 'trans_y', 'trans_z', 
+                          'rot_x', 'rot_y', 'rot_z']] # can add 'global_signal' also ,,
+                          # 
+     # change NaN of FD to zero
+    finalConf = np.array(finalConf.fillna(0.0))
+    return finalConf
+
+
+# %% functional files
+
+func_template = '/gpfs/gibbs/pi/levy_ifat/Or/kpe/fmriprep/sub-{sub}/ses-{session}/func/sub-{sub}_ses-{session}_task-rest_space-MNI152NLin2009cAsym_desc-preproc_bold.nii.gz'
+confound_template = '/gpfs/gibbs/pi/levy_ifat/Or/kpe/fmriprep/sub-{sub}/ses-{session}/func/sub-{sub}_ses-{session}_task-rest_desc-confounds_regressors.tsv'
+
+# get subject list
+# from bids import BIDSLayout
+# folder = '/gpfs/gibbs/pi/levy_ifat/Or/kpe/fmriprep/'
+# layout = BIDSLayout(folder, validate=False)
+# subject_list = layout.get_subjects()
+# %%
+# create a mean mask of all subjects
+# load mask of brain
+
+
+brainmasks = glob.glob('/gpfs/gibbs/pi/levy_ifat/Or/kpe/fmriprep/sub-*/ses-%s/func/sub-*_ses-%s_task-rest_space-MNI152NLin2009cAsym_desc-brain_mask.nii.gz' %(ses,ses))
+print(brainmasks)
+# %matplotlib inline
+#for mask in brainmasks:
+ #   nilearn.plotting.plot_roi(mask)
+
+mean_mask = nilearn.image.mean_img(brainmasks)
+#nilearn.plotting.plot_stat_map(mean_mask)
+group_mask = nilearn.image.math_img("a>=0.98", a=mean_mask)
+#nilearn.plotting.plot_roi(group_mask)
+
+# %% fetch atlas
+maps_img = '/gpfs/gibbs/pi/levy_ifat/Or/DiFuMo_atlas/256/maps.nii.gz'
+labels = pd.read_csv('/gpfs/gibbs/pi/levy_ifat/Or/DiFuMo_atlas/256/labels_256_dictionary.csv')
+#coords = nilearn.plotting.find_parcellation_cut_coords(labels_img=maps_img)
+coords = nilearn.plotting.find_probabilistic_atlas_cut_coords(maps_img)
+# generate time series
+#
+mask_params = { 'mask_img': group_mask,
+               'detrend': True, 'standardize': True,
+               'high_pass': 0.01, 'low_pass': 0.1, 't_r': 1,
+               'smoothing_fwhm': 6.,
+                'verbose': 5}
+
+masker = NiftiMapsMasker(maps_img=maps_img, **mask_params)
+
+# %%
+subject_list = ['1587']
+
+# %% Generate npy files of timeseries for each subject per session
+# we will use it later on, stratify to scripts etc.
+# build a specific folder
+try:
+    os.makedirs(output_dir)
+except:
+    print('Folder already exist')
+
+subject_ts = []
+for sub in subject_list:
+    print(f' Analysing subject {sub}')
+    subject = sub
+    func = func_template.format(sub=subject, session=ses)
+    confound = confound_template.format(sub=subject, session=ses)
+    try:
+        signals = masker.fit_transform(imgs=func, confounds=removeVars(confound))
+        save = np.save(output_dir + 'sub-' + subject + '_ses-' + ses, signals)
+        subject_ts.append(signals)
+    except:
+        print(f'Subject {sub} has no data')
+
+
+
+
+
+# %%