Make p-value adjustment optional; bug fixes for negative valued clusters

Author: davidmwatson

defineROICluster_Freesurfer.py

@@ -15,7 +15,8 @@
     Surface hemisphere: 'lh' or 'rh'
 
 -i|--input
-    Path to input statistical surface image (probably a -log10(pstat))
+    Path to input statistical surface image (probably a -log10(pstat) or
+    zstat).
 
 -n|--cluster-size
     Desired size of cluster, in units of mm^2 or number of vertices depending
@@ -39,7 +40,7 @@
     ROIs - the number must match the number of seed coordinates.
 
 --surf (optional)
-    Which surface to use (default = white)
+    Which surface to use (default = white). Note: do not use inflated surface.
 
 -l|--initial-label (optional)
     Path to surface label file. If provided, will apply this label to the data
@@ -67,11 +68,16 @@
 
 --sign (optional)
     Tail of statistical map to use. Can be abs, pos (default), or neg. Note
-    that thresholds will be adjusted if a one-tailed pos/neg test is used, so
-    that values will correspond to one-tailed p-values. For example, if sign is
-    'pos' and the threshold is reported as 3.0 (one-tailed p = .001), the
-    statistical map would actually be thresholded at 3.0 - log(2) ~= 2.7
-    (two-tailed p = .002).
+    that for a negative tail, all thresholds and bounds should still be
+    specified as positive values (they will be sign-flipped internally).
+
+--adjust-pvals (optional)
+    If specified, adjust threshold for one-tailed tests (ignored if --sign is
+    'abs'). Assumes that input is a -log10(pstat) image. Values will be
+    adjusted by -log10(2) (equivalent to dividing the p-value by 2). For
+    example, if sign is 'pos' and the threshold is reported as 3.0 (one-tailed
+    p = .001), the input image would actually be thresholded at
+    3.0 - log(2) ~= 2.7 (two-tailed p = .002).
 
 --subjects-dir (optional)
     Directory containing Freesurfer subject directories. If omitted, defaults
@@ -130,7 +136,7 @@ def my_subprocess_run(*args, **kwargs):
     return rtn
 
 
-def RAS_to_vertex(subjects_dir, subject, hemi, vtx_coord):
+def RAS_to_vertex(subjects_dir, subject, hemi, surf, vtx_coord):
     """
     Get index of closest vertex to surfaceRAS co-ordinate.
 
@@ -142,23 +148,24 @@ def RAS_to_vertex(subjects_dir, subject, hemi, vtx_coord):
         Subject ID
     hemi : str
         lh or rh
+    surf : str
+        Surface to use (e.g., white)
     vtx_coord : array-like
         [x, y, z] surface RAS co-ordindate of vertex
 
-
     Returns
     -------
     vtx : int
         Index of vertex
     """
-    geom_file = os.path.join(subjects_dir, subject, 'surf', hemi + '.white')
+    geom_file = os.path.join(subjects_dir, subject, 'surf', f'{hemi}.{surf}')
     surf_coords = nib.freesurfer.read_geometry(geom_file)[0]
     vtx_coord = np.asarray(vtx_coord, dtype=float)
     return np.linalg.norm(surf_coords - vtx_coord, axis=1).argmin()
 
 
 def mri_surfcluster(thr, subjects_dir, subject, hemi, surf, infile,
-                    clusterdir, sign, clabel=None, save_labels=False):
+                    clusterdir, sign, clabel, save_labels, adjust_pvals):
     """
     Run Freesurfer mri_surfcluster command via subprocess
 
@@ -181,10 +188,13 @@ def mri_surfcluster(thr, subjects_dir, subject, hemi, surf, infile,
         Path to directory to save outputs. Must already exist.
     sign : str
         abs, pos, or neg
-    clabel : str (optional)
-        Path to label file to apply as pre-mask
-    save_labels : bool (optional)
+    clabel : str or None
+        Optional path to label file to apply as pre-mask
+    save_labels : bool
         If True, also save clusters out as label files. Default is False.
+    adjust_pvals : bool
+        If True, assume threshold is for a -log10(pstat), and adjust for
+        one-tailed test.
 
     Outputs
     -------
@@ -213,6 +223,8 @@ def mri_surfcluster(thr, subjects_dir, subject, hemi, surf, infile,
         cmd.extend(['--clabel', clabel])
     if save_labels:
         cmd.extend(['--olab', os.path.join(clusterdir, 'cluster')])
+    if not adjust_pvals:
+        cmd.extend(['--no-adjust'])
 
     # Run command
     my_subprocess_run(cmd)
@@ -290,24 +302,26 @@ def get_cluster_size(clusterdir, clusterID, size_metric):
     raise Exception('Cluster ID not found')
 
 
-def cost_func(thr, subjects_dir, subject, hemi, surf, infile, clusterdir,
-              sign, seed_vertex, target_size, size_metric, clabel=None):
+def cost_func(thr, seed_vertex, target_size, size_metric,
+              mri_surfcluster_kwargs):
     """
     Cost function to be optimised. Applies clustering, then returns error
     between actual and target cluster size.
 
     Arguments
     ---------
-    thr, subjects_dir, subject, hemi, surf, clusterdir, surf, clabel
-        All as per mri_surfcluster function. Threshold may be supplied as
-        single-item array - this is for compatibility with minimize function.
+    thr : float
+        Threshold to apply. May be supplied as single-item array - this is for
+        compatibility with minimize function.
     seed_vertex : int
         Index of seed vertex
     target_size : float or int
         Target size of cluster; either area (in mm^2) or num vertices
         depending on size metric.
     size_metric : str
         Size metric: 'area' or 'nVtxs'
+    mri_surfcluster_kwargs : dict
+        Dictionary of further keyword arguments to be passed to mri_surfcluster
 
     Returns
     -------
@@ -322,10 +336,10 @@ def cost_func(thr, subjects_dir, subject, hemi, surf, infile, clusterdir,
             raise ValueError('thr must be float or single-item array')
 
     # Cluster
-    mri_surfcluster(thr, subjects_dir, subject, hemi, surf, infile, clusterdir,
-                    sign, clabel=clabel)
+    mri_surfcluster(thr, **mri_surfcluster_kwargs)
 
     # Find size of cluster
+    clusterdir = mri_surfcluster_kwargs['clusterdir']
     clusterID = get_cluster_at_vertex(clusterdir, seed_vertex)
     if clusterID == 0:
         cluster_size = 0
@@ -375,6 +389,8 @@ class CustomFormatter(argparse.RawTextHelpFormatter,
                     help='Upper-bound threshold; set to inf or nan to disable')
 parser.add_argument('--sign', choices=['abs', 'pos', 'neg'], default='pos',
                     help='Tail of statistical map to use')
+parser.add_argument('--adjust-pvals', action='store_true',
+                    help='Adjust threshold for one-tailed p-values')
 parser.add_argument('--subjects-dir', help='Freesurfer subjects directory')
 
 # If no arguments given, print help and exit
@@ -397,6 +413,7 @@ class CustomFormatter(argparse.RawTextHelpFormatter,
 lower_bound = args.lower_bound
 upper_bound = args.upper_bound
 sign = args.sign
+adjust_pvals = args.adjust_pvals
 subjects_dir = args.subjects_dir
 
 # Allocate default subjects dir if necessary
@@ -433,7 +450,7 @@ class CustomFormatter(argparse.RawTextHelpFormatter,
     # If seed in scannerRAS coords, convert to vertex index. If already index,
     # extract value and convert float -> int
     if len(seed_vertex) == 3:
-        seed_vertex = RAS_to_vertex(subjects_dir, subject, hemi, seed_vertex)
+        seed_vertex = RAS_to_vertex(subjects_dir, subject, hemi, surf, seed_vertex)
     elif len(seed_vertex) == 1:
         seed_vertex = int(seed_vertex[0])
     else:
@@ -446,11 +463,30 @@ class CustomFormatter(argparse.RawTextHelpFormatter,
     # Check lower bound appropriate for seed
     if lower_bound is not None:
         seed_val = surf_data[seed_vertex]
-        if sign in ['pos','neg']:
+
+        if sign == 'abs':
+            seed_val = abs(seed_val)
+        elif sign == 'neg':
+            seed_val = -1 * seed_val
+
+        if sign in ['pos','neg'] and adjust_pvals:
             seed_val -= np.log10(2)
+
         if seed_val < lower_bound:
             raise ValueError('Lower bound cannot exceed value at seed voxel')
 
+    # Build dict of kwargs for mri_surfcluster
+    mri_surfcluster_kwargs = {'subjects_dir':subjects_dir,
+                              'subject':subject,
+                              'hemi':hemi,
+                              'surf':surf,
+                              'infile':infile,
+                              'clusterdir':tmpdir.name,
+                              'sign':sign,
+                              'clabel':initial_label,
+                              'save_labels':False,
+                              'adjust_pvals':adjust_pvals}
+
     # Loop starting parameters
     all_thr = []
     all_errs = []
@@ -459,17 +495,23 @@ class CustomFormatter(argparse.RawTextHelpFormatter,
         # Otherwise, just convert to float.
         if x0.endswith('%'):
             prop = float(x0.rstrip('%')) / 100
-            max_ = surf_data[seed_vertex]
+
             min_ = lower_bound if lower_bound is not None else 0
+            max_ = surf_data[seed_vertex]
+            if sign == 'abs':
+                max_ = abs(max_)
+            elif sign == 'neg':
+                max_ = -1 * max_
+
             x0 = min_ + prop * (max_ - min_)
-            if sign in ['pos','neg']:
+            if sign in ['pos','neg'] and adjust_pvals:
                 x0 += np.log10(2)
         else:
             x0 = float(x0)
 
+
         # Run optimisation
-        args = (subjects_dir, subject, hemi, surf, infile, tmpdir.name, sign,
-                seed_vertex, target_size, size_metric, initial_label)
+        args = (seed_vertex, target_size, size_metric, mri_surfcluster_kwargs)
         optimres = minimize(cost_func, x0=[x0], args=args, method='Nelder-Mead')
 
         # Extract optimised threshold and error values, append to array
@@ -489,8 +531,8 @@ class CustomFormatter(argparse.RawTextHelpFormatter,
         opt_thr = upper_bound
 
     # Run clustering once more with selected thresh
-    mri_surfcluster(opt_thr, subjects_dir, subject, hemi, surf, infile,
-                    tmpdir.name, sign, clabel=initial_label, save_labels=True)
+    mri_surfcluster_kwargs['save_labels'] = True
+    mri_surfcluster(opt_thr, **mri_surfcluster_kwargs)
     clusterID = get_cluster_at_vertex(tmpdir.name, seed_vertex)
     cluster_size = get_cluster_size(tmpdir.name, clusterID, size_metric)
 
@@ -505,4 +547,3 @@ class CustomFormatter(argparse.RawTextHelpFormatter,
 
 # Cleanup tmpdir
 tmpdir.cleanup()
-