Merge pull request #181 from CoBrALab/updates

Updates

Dockerfile

@@ -50,7 +50,7 @@ RUN mkdir -p /opt/ANTs/build && git clone https://github.com/ANTsX/ANTs.git /opt
     && cd /opt/ANTs/src \
     && git checkout 1759e5e23772e114a490cfa33a5764b400307b9d \
     && cd /opt/ANTs/build \
-    && cmake -GNinja -DITK_BUILD_MINC_SUPPORT=ON ../src \
+    && cmake -GNinja -DITK_BUILD_MINC_SUPPORT=ON -DBUILD_TESTING=OFF ../src \
     && cmake --build . \
     && cd ANTS-build \
     && cmake --install .

README.md

@@ -588,16 +588,16 @@ rabies -p MultiProc preprocess test_dataset/ preprocess_outputs/ --TR 1.0s --no_
 First, this will run the minimal preprocessing step on the test dataset and store outputs into preprocess_outputs/ folder. The option -p MultiProc specifies to run the pipeline in parallel according to available local threads.
 <br/>
 
-**confound_regression**
+**confound_correction**
 ```sh
-rabies -p MultiProc confound_regression preprocess_outputs/ confound_regression_outputs/ --TR 1.0s --commonspace_bold --smoothing_filter 0.3 --conf_list WM_signal CSF_signal vascular_signal mot_6
+rabies -p MultiProc confound_correction preprocess_outputs/ confound_correction_outputs/ --TR 1.0s --commonspace_bold --smoothing_filter 0.3 --conf_list WM_signal CSF_signal vascular_signal mot_6
 ```
-Next, to conduct the modeling and regression of confounding sources, the confound_regression step can be run with custom options for denoising. In this case, we apply a highpass filtering at 0.01Hz, together with the voxelwise regression of the 6 rigid realignment parameters and the mean WM,CSF and vascular signal which are derived from masks provided along with the anatomical template. Finally, a smoothing filter 0.3mm is applied. We are running this on the commonspace outputs from preprocess (--commonspace_bold), since we will run analysis in commonspace in the next step.
+Next, to conduct the modeling and regression of confounding sources, the confound_correction step can be run with custom options for denoising. In this case, we apply a highpass filtering at 0.01Hz, together with the voxelwise regression of the 6 rigid realignment parameters and the mean WM,CSF and vascular signal which are derived from masks provided along with the anatomical template. Finally, a smoothing filter 0.3mm is applied. We are running this on the commonspace outputs from preprocess (--commonspace_bold), since we will run analysis in commonspace in the next step.
 <br/>
 
 **analysis**
 ```sh
-rabies -p MultiProc analysis confound_regression_outputs analysis_outputs/ --TR 1.0s --group_ICA --DR_ICA
+rabies -p MultiProc analysis confound_correction_outputs analysis_outputs/ --TR 1.0s --group_ICA --DR_ICA
 ```
 Finally, RABIES has a few standard analysis options provided, which are specified in the Analysis documentation. In this example, we are going to run group independent component analysis (--group_ICA), using FSL's MELODIC function, followed by a dual regression (--DR_ICA) to back propagate the group components onto individual subjects.
 
@@ -627,12 +627,12 @@ singularity run -B $PWD/test_dataset:/test_dataset:ro \
 ```
 <br/>
 
-**confound_regression**
+**confound_correction**
 ```sh
 singularity run -B $PWD/test_dataset:/test_dataset:ro \
 -B $PWD/preprocess_outputs:/preprocess_outputs/ \
--B $PWD/confound_regression_outputs:/confound_regression_outputs/ \
-/path_to_singularity_image/rabies.sif -p MultiProc confound_regression /preprocess_outputs/ /confound_regression_outputs/ --TR 1.0s --highpass 0.01 --commonspace_bold --smoothing_filter 0.3 --conf_list WM_signal CSF_signal vascular_signal mot_6
+-B $PWD/confound_correction_outputs:/confound_correction_outputs/ \
+/path_to_singularity_image/rabies.sif -p MultiProc confound_correction /preprocess_outputs/ /confound_correction_outputs/ --TR 1.0s --highpass 0.01 --commonspace_bold --smoothing_filter 0.3 --conf_list WM_signal CSF_signal vascular_signal mot_6
 ```
 Note here that the path to the dataset is still linked to the container with -B, even though it is not explicitely part of the inputs in the confound regression call. This is necessary since the paths used in the preprocess steps are still accessed in the background, and there will be an error if the paths are not kept consistent across processing steps.
 <br/>
@@ -641,9 +641,9 @@ Note here that the path to the dataset is still linked to the container with -B,
 ```sh
 singularity run -B $PWD/test_dataset:/test_dataset:ro \
 -B $PWD/preprocess_outputs:/preprocess_outputs/ \
--B $PWD/confound_regression_outputs:/confound_regression_outputs/ \
+-B $PWD/confound_correction_outputs:/confound_correction_outputs/ \
 -B $PWD/analysis_outputs:/analysis_outputs/ \
-/path_to_singularity_image/rabies.sif -p MultiProc analysis /confound_regression_outputs /analysis_outputs/ --TR 1.0s --group_ICA --DR_ICA
+/path_to_singularity_image/rabies.sif -p MultiProc analysis /confound_correction_outputs /analysis_outputs/ --TR 1.0s --group_ICA --DR_ICA
 ```
 <br/>
 
@@ -733,8 +733,8 @@ The following image presents an example of the overlap for the EPI2Anat registra
 <details><summary><b>Click to expand</b></summary>
 <p>
 
-Important outputs from confound regression will be found in the confound_regression_datasink present in the provided output folder:
-- **confound_regression_datasink**: Includes outputs specific to the anatomical preprocessing workflow
+Important outputs from confound regression will be found in the confound_correction_datasink present in the provided output folder:
+- **confound_correction_datasink**: Includes outputs specific to the anatomical preprocessing workflow
     - cleaned_timeseries: Resulting timeseries after the application of confound regression
     - VE_file: .pkl file which contains a dictionary vectors, where each vector corresponds to the voxelwise the variance explained (VE) from each regressor in the regression model
     - aroma_out: if --run_aroma is selected, the outputs from running ICA-AROMA will be saved, which includes the MELODIC ICA outputs and the component classification results

rabies/__version__.py

@@ -3,6 +3,6 @@
 # 88YbdP88   8P   88"""   dP__Yb  Yb      88"Yb   dP__Yb  Yb  "88 88""
 # 88 YY 88  dP    88     dP""""Yb  YboodP 88  Yb dP""""Yb  YboodP 888888
 
-VERSION = (0, 3, 2)
+VERSION = (0, 3, 3)
 
 __version__ = '.'.join(map(str, VERSION))

rabies/analysis_pkg/analysis_functions.py

@@ -1,6 +1,6 @@
 import numpy as np
 import nibabel as nb
-
+from .analysis_math import vcorrcoef,closed_form
 
 def seed_based_FC(bold_file, brain_mask, seed_dict, seed_name):
     import os
@@ -57,15 +57,6 @@ def seed_corr(bold_file, brain_mask, seed):
     return corrs
 
 
-def vcorrcoef(X, y):  # return a correlation between each row of X with y
-    Xm = np.reshape(np.mean(X, axis=1), (X.shape[0], 1))
-    ym = np.mean(y)
-    r_num = np.sum((X-Xm)*(y-ym), axis=1)
-    r_den = np.sqrt(np.sum((X-Xm)**2, axis=1)*np.sum((y-ym)**2))
-    r = r_num/r_den
-    return r
-
-
 def get_CAPs(data, volumes, n_clusters):
     from sklearn.cluster import KMeans
     kmeans = KMeans(n_clusters=n_clusters, n_init=10, max_iter=300)
@@ -209,7 +200,7 @@ def plot_matrix(filename, corr_matrix):
 '''
 
 
-def run_group_ICA(bold_file_list, mask_file, dim):
+def run_group_ICA(bold_file_list, mask_file, dim, random_seed):
     import os
     import pandas as pd
 
@@ -222,7 +213,7 @@ def run_group_ICA(bold_file_list, mask_file, dim):
 
     from rabies.preprocess_pkg.utils import run_command
     out_dir = os.path.abspath('group_melodic.ica')
-    command = f'melodic -i {file_path} -m {mask_file} -o {out_dir} -d {dim} --report --seed=1'
+    command = f'melodic -i {file_path} -m {mask_file} -o {out_dir} -d {dim} --report --seed={str(random_seed)}'
     rc = run_command(command)
     IC_file = out_dir+'/melodic_IC.nii.gz'
     return out_dir, IC_file
@@ -307,21 +298,6 @@ def run_DR_ICA(bold_file, mask_file, IC_file):
     return DR_maps_filename, dual_regression_timecourse_csv
 
 
-'''
-LINEAR REGRESSION --- CLOSED-FORM SOLUTION
-'''
-
-
-def closed_form(X, Y, intercept=False):  # functions that computes the Least Squares Estimates
-    if intercept:
-        X = np.concatenate((X, np.ones([X.shape[0], 1])), axis=1)
-    return np.linalg.inv(X.transpose().dot(X)).dot(X.transpose()).dot(Y)
-
-
-def mse(X, Y, w):  # function that computes the Mean Square Error (MSE)
-    return np.mean((Y-np.matmul(X, w))**2)
-
-
 def dual_regression(all_IC_vectors, timeseries):
     ### compute dual regression
     ### Here, we adopt an approach where the algorithm should explain the data

rabies/analysis_pkg/analysis_math.py

@@ -0,0 +1,101 @@
+import numpy as np
+
+def vcorrcoef(X, y):  # return a correlation between each row of X with y
+    Xm = np.reshape(np.mean(X, axis=1), (X.shape[0], 1))
+    ym = np.mean(y)
+    r_num = np.sum((X-Xm)*(y-ym), axis=1)
+    r_den = np.sqrt(np.sum((X-Xm)**2, axis=1)*np.sum((y-ym)**2))
+    r = r_num/r_den
+    return r
+
+
+def elementwise_corrcoef(X, Y):
+    # X and Y are each of shape num_observations X num_element
+    # computes the correlation between each element of X and Y
+    Xm = X.mean(axis=0)
+    Ym = Y.mean(axis=0)
+    r_num = np.sum((X-Xm)*(Y-Ym), axis=0)
+
+    r_den = np.sqrt(np.sum((X-Xm)**2, axis=0)*np.sum((Y-Ym)**2, axis=0))
+    r = r_num/r_den
+    return r
+
+
+def elementwise_spearman(X,Y):    
+    order = X.argsort(axis=0)
+    X_ranks = order.argsort(axis=0)
+    order = Y.argsort(axis=0)
+    Y_ranks = order.argsort(axis=0)
+    return elementwise_corrcoef(X_ranks, Y_ranks)
+
+
+def dice_coefficient(mask1,mask2):
+    dice = np.sum(mask1*mask2)*2.0 / (np.sum(mask1) + np.sum(mask2))
+    return dice
+
+
+'''
+LINEAR REGRESSION --- CLOSED-FORM SOLUTION
+'''
+
+
+def closed_form(X, Y, intercept=False):  # functions that computes the Least Squares Estimates
+    if intercept:
+        X = np.concatenate((X, np.ones([X.shape[0], 1])), axis=1)
+    return np.linalg.inv(X.transpose().dot(X)).dot(X.transpose()).dot(Y)
+
+
+def mse(X, Y, w):  # function that computes the Mean Square Error (MSE)
+    return np.mean((Y-np.matmul(X, w))**2)
+
+
+
+'''
+
+def closed_form_3d(X,Y):
+    return np.matmul(np.matmul(np.linalg.inv(np.matmul(X.transpose(0,2,1),X)),X.transpose(0,2,1)),Y)
+
+def lme_stats_3d(X,Y):
+    #add an intercept
+    X=np.concatenate((X,np.ones((X.shape[0],X.shape[1],1))),axis=2)
+    [num_comparisons,num_observations,num_predictors] = X.shape
+    [num_comparisons,num_observations,num_features] = Y.shape
+
+    w=closed_form_3d(X,Y)
+
+    residuals = Y-np.matmul(X, w)
+    MSE = (((residuals)**2).sum(axis=1)/(num_observations-num_predictors))
+
+
+    var_b = np.expand_dims(MSE, axis=1)*np.expand_dims(np.linalg.inv(np.matmul(X.transpose(0,2,1),X)).diagonal(axis1=1,axis2=2), axis=2)
+    sd_b = np.sqrt(var_b) # standard error on the Betas
+    ts_b = w/sd_b # calculate t-values for the Betas
+    p_values =[2*(1-stats.t.cdf(np.abs(ts_b[:,i,:]),(num_observations-num_predictors))) for i in range(ts_b.shape[1])] # calculate a p-value map for each predictor
+
+    return ts_b,p_values,w,residuals
+
+non_nan_idx = (np.isnan(voxelwise_array).sum(axis=(0,1))==0)
+
+# take out the voxels which have null values
+X=voxelwise_array[:,:6,non_nan_idx].transpose(2,0,1)
+Y=voxelwise_array[:,6:,non_nan_idx].transpose(2,0,1)
+
+
+ts_b,p_values,w,residuals = lme_stats_3d(X,Y)
+
+x_name=['Somatomotor','Dorsal Comp','DMN', 'Prior Modeling 1', 'Prior Modeling 2', 'Prior Modeling 3']
+y_name=['group','temporal_std','VE_spatial','GS_corr','DVARS_corr','FD_corr']
+
+fig,axes = plt.subplots(nrows=len(x_name), ncols=len(y_name),figsize=(12*len(y_name),3*len(x_name)))
+
+
+for i,x_label in zip(range(len(x_name)),x_name):
+    for j,y_label in zip(range(len(y_name)),y_name):
+        ax=axes[i,j]
+
+        stat_map=np.zeros(voxelwise_array.shape[2])
+        stat_map[non_nan_idx]=ts_b[:,j,i]
+
+        ax.set_title('T-value of {} on {}'.format(y_label,x_label), fontsize=15)
+        plot_stat_map(analysis_functions.recover_3D(mask_file,stat_map),bg_img='DSURQE.nii.gz', axes=ax, cut_coords=(0,1,2,3,4,5), display_mode='y')
+'''

rabies/analysis_pkg/analysis_wf.py

@@ -61,7 +61,7 @@ def init_analysis_wf(opts, commonspace_cr=False, seed_list=[], name="analysis_wf
 
     include_group_ICA = opts.group_ICA
 
-    if opts.DR_ICA and not opts.data_diagnosis:
+    if opts.DR_ICA or opts.data_diagnosis:
         if not commonspace_cr:
             raise ValueError(
                 'Outputs from confound regression must be in commonspace to run dual regression. Try running confound regression again with --commonspace_analysis.')
@@ -91,11 +91,12 @@ def init_analysis_wf(opts, commonspace_cr=False, seed_list=[], name="analysis_wf
         if not commonspace_cr:
             raise ValueError(
                 'Outputs from confound regression must be in commonspace to run group-ICA. Try running confound regression again with --nativespace_analysis.')
-        group_ICA = pe.Node(Function(input_names=['bold_file_list', 'mask_file', 'dim'],
+        group_ICA = pe.Node(Function(input_names=['bold_file_list', 'mask_file', 'dim', 'random_seed'],
                                      output_names=['out_dir', 'IC_file'],
                                      function=run_group_ICA),
                             name='group_ICA', mem_gb=1)
         group_ICA.inputs.dim = opts.dim
+        group_ICA.inputs.random_seed = opts.melodic_random_seed
 
         workflow.connect([
             (group_inputnode, group_ICA, [
@@ -108,16 +109,11 @@ def init_analysis_wf(opts, commonspace_cr=False, seed_list=[], name="analysis_wf
                 ]),
             ])
 
-    if opts.dual_ICA>0 and not opts.data_diagnosis:
+    if opts.dual_ICA>0:
         if not commonspace_cr:
             raise ValueError(
                 'Outputs from confound regression must be in commonspace to run group-ICA. Try running confound regression again with --nativespace_analysis.')
         from .analysis_functions import dual_ICA_wrapper
-        from .data_diagnosis import resample_IC_file
-        resample_IC = pe.Node(Function(input_names=['in_file', 'ref_file'],
-                                     output_names=['out_file'],
-                                     function=resample_IC_file),
-                            name='resample_IC', mem_gb=1)
 
         dual_ICA_node = pe.Node(dual_ICA_wrapper(),
                             name='dual_ICA_node', mem_gb=1)