diff --git a/pynets/cli/pynets_collect.py b/pynets/cli/pynets_collect.py index 68a916c7..152d9f71 100644 --- a/pynets/cli/pynets_collect.py +++ b/pynets/cli/pynets_collect.py @@ -420,7 +420,8 @@ def recover_missing(bad_col, bad_cols_dict, rerun_dict, modality, else: # from pynets.stats.netstats import collect_pandas_df_make # collect_pandas_df_make(glob.glob(f"{working_path}/{sub}/{ses}/" - # f"{modality}/{atlas}/topology/*_neat.csv"), + # f"{modality}/{atlas}/ + # topology/*_neat.csv"), # f"{sub}_{ses}", None, False) rerun_dict[sub][ses][modality][atlas].append(bad_col) continue @@ -461,7 +462,8 @@ def recover_missing(bad_col, bad_cols_dict, rerun_dict, modality, # collect_pandas_df_make # collect_pandas_df_make( # glob.glob(f"{working_path}/{sub}/{ses}/" - # f"{modality}/{atlas}/topology/*_neat.csv"), + # f"{modality}/{atlas}/topology/ + # *_neat.csv"), # f"{sub}_{ses}", None, False) continue del df_tmp @@ -573,8 +575,8 @@ def is_non_zero_file(fpath): set( [ os.path.basename(str(Path(i).parent.parent)) - for i in glob.glob(f"{working_path}/{sub}/*/{modality}/*/topology/*", - recursive=True) + for i in glob.glob(f"{working_path}/{sub}/*/{modality}/*/" + f"topology/*", recursive=True) ] ) ) @@ -894,7 +896,7 @@ def main(): import sys import glob from pynets.cli.pynets_collect import build_collect_workflow - from types import SimpleNamespace + # from types import SimpleNamespace from pathlib import Path try: @@ -910,28 +912,28 @@ def main(): " flag.\n") sys.exit() - # args = get_parser().parse_args() - args_dict_all = {} - args_dict_all['plug'] = 'MultiProc' - args_dict_all['v'] = False - args_dict_all['pm'] = '24,57' - #args_dict_all['basedir'] = '/working/tuning_set/outputs_shaeffer/pynets' - #args_dict_all['basedir'] = '/scratch/04171/dpisner/HNU/HNU_outs/triple/pynets' - args_dict_all['basedir'] = '/scratch/04171/dpisner/HNU/HNU_outs/visual/pynets' - args_dict_all['work'] = '/tmp/work/dwi' - args_dict_all['modality'] = 'dwi' - args_dict_all['dc'] = ['diversity_coefficient', - 'participation_coefficient', - 'average_local_efficiency', - 'average_clustering', - 'average_local_clustering_nodewise', - 'average_local_efficiency_nodewise', - 'degree_centrality', - 'weighted_transitivity', - # "_minlength-0", - "_minlength-20", "_minlength-30", "variance", - "res-1000"] - args = SimpleNamespace(**args_dict_all) + args = get_parser().parse_args() + # args_dict_all = {} + # args_dict_all['plug'] = 'MultiProc' + # args_dict_all['v'] = False + # args_dict_all['pm'] = '24,57' + # #args_dict_all['basedir'] = '/working/tuning_set/outputs_shaeffer/pynets' + # #args_dict_all['basedir'] = '/scratch/04171/dpisner/HNU/HNU_outs/triple/pynets' + # args_dict_all['basedir'] = '/scratch/04171/dpisner/HNU/HNU_outs/visual/pynets' + # args_dict_all['work'] = '/tmp/work/dwi' + # args_dict_all['modality'] = 'dwi' + # args_dict_all['dc'] = ['diversity_coefficient', + # 'participation_coefficient', + # 'average_local_efficiency', + # 'average_clustering', + # 'average_local_clustering_nodewise', + # 'average_local_efficiency_nodewise', + # 'degree_centrality', + # 'weighted_transitivity', + # # "_minlength-0", + # "_minlength-20", "_minlength-30", "variance", + # "res-1000"] + # args = SimpleNamespace(**args_dict_all) from multiprocessing import set_start_method, Process, Manager @@ -954,9 +956,13 @@ def main(): gc.collect() mgr.shutdown() - working_path = args_dict_all['basedir'] - modality = args_dict_all['modality'] - drop_cols = args_dict_all['dc'] + # working_path = args_dict_all['basedir'] + # modality = args_dict_all['modality'] + # drop_cols = args_dict_all['dc'] + working_path = args.basedir + modality = args.modality + drop_cols = args.dc + all_files = glob.glob( f"{str(Path(working_path))}/{modality}_group_topology_auc/*.csv" ) diff --git a/pynets/cli/pynets_run.py b/pynets/cli/pynets_run.py index 81509253..0c6a2fc1 100755 --- a/pynets/cli/pynets_run.py +++ b/pynets/cli/pynets_run.py @@ -174,7 +174,7 @@ def get_parser(): metavar="Path to custom parcellation file", default=None, nargs="+", - help="(Hyperparameter): Optionally specify a path to a " + help="(metaparameter): Optionally specify a path to a " "parcellation/atlas Nifti1Image file in MNI space. Labels should" "be spatially distinct across hemispheres and ordered with " "consecutive integers with a value of 0 as the background label." @@ -182,7 +182,7 @@ def get_parser(): "separate them by space.\n", ) - # Modality-pervasive hyperparameters + # Modality-pervasive metaparameters parser.add_argument( "-mod", metavar="Connectivity estimation/reconstruction method", @@ -201,15 +201,14 @@ def get_parser(): "csd", "sfm", ], - help="(Hyperparameter): Specify connectivity estimation model. " + help="(metaparameter): Specify connectivity estimation model. " "For fMRI, possible models include: corr for correlation, " "cov for covariance, sps for precision covariance, partcorr for " - "partial correlation. sps type is used by default. If skgmm is " + "partial correlation. If skgmm is " "installed (https://github.com/skggm/skggm), then " "QuicGraphicalLasso, QuicGraphicalLassoCV, " - "QuicGraphicalLassoEBIC, and AdaptiveQuicGraphicalLasso. Default " - "is partcorr for fMRI. For dMRI, current models include csa, " - "csd, and sfm.\n", + "QuicGraphicalLassoEBIC, and AdaptiveQuicGraphicalLasso. For " + "dMRI, current models include csa, csd, and sfm.\n", ) parser.add_argument( "-a", @@ -251,7 +250,7 @@ def get_parser(): 'sub-colin27_label-L2018_desc-scale4_atlas', 'sub-colin27_label-L2018_desc-scale5_atlas' ], - help="(Hyperparameter): Specify an atlas parcellation from nilearn or " + help="(metaparameter): Specify an atlas parcellation from nilearn or " "local libraries. If you wish to iterate your pynets run over " "multiple atlases, separate them by space. Available nilearn " "atlases are:" @@ -280,7 +279,7 @@ def get_parser(): metavar="Spherical centroid node size", default=4, nargs="+", - help="(Hyperparameter): Optionally specify coordinate-based node " + help="(metaparameter): Optionally specify coordinate-based node " "radius size(s). Default is 4 " "mm for fMRI and 8mm for dMRI. If you wish to iterate the pipeline " "across multiple node sizes, separate the list " @@ -298,29 +297,29 @@ def get_parser(): "-min_thr", metavar="Multi-thresholding minimum threshold", default=None, - help="(Hyperparameter): Minimum threshold for multi-thresholding.\n", + help="(metaparameter): Minimum threshold for multi-thresholding.\n", ) parser.add_argument( "-max_thr", metavar="Multi-thresholding maximum threshold", default=None, - help="(Hyperparameter): Maximum threshold for multi-thresholding.", + help="(metaparameter): Maximum threshold for multi-thresholding.", ) parser.add_argument( "-step_thr", metavar="Multi-thresholding step size", default=None, - help="(Hyperparameter): Threshold step value for multi-thresholding. " + help="(metaparameter): Threshold step value for multi-thresholding. " "Default is 0.01.\n", ) - # fMRI hyperparameters + # fMRI metaparameters parser.add_argument( "-sm", metavar="Smoothing value (mm fwhm)", default=0, nargs="+", - help="(Hyperparameter): Optionally specify smoothing width(s). " + help="(metaparameter): Optionally specify smoothing width(s). " "Default is 0 / no smoothing. If you wish to iterate the pipeline" " across multiple smoothing separate the list " "by space (e.g. 2 4 6).\n", @@ -330,7 +329,7 @@ def get_parser(): metavar="High-pass filter (Hz)", default=None, nargs="+", - help="(Hyperparameter): Optionally specify high-pass filter values " + help="(metaparameter): Optionally specify high-pass filter values " "to apply to node-extracted time-series for fMRI. " "Default is None. If you wish to iterate the pipeline across " "multiple high-pass filter thresholds, values, " @@ -361,7 +360,7 @@ def get_parser(): metavar="Number of k clusters", default=None, nargs="+", - help="(Hyperparameter): Specify a number of clusters to produce. " + help="(metaparameter): Specify a number of clusters to produce. " "If you wish to iterate the pipeline across multiple values of k," " separate the list by space (e.g. 100 150 200).\n", ) @@ -372,7 +371,7 @@ def get_parser(): nargs="+", choices=["ward", "rena", "kmeans", "complete", "average", "single", "ncut"], - help="(Hyperparameter): Specify the types of clustering to use. " + help="(metaparameter): Specify the types of clustering to use. " "Recommended options are: ward, rena, kmeans, or ncut. Note that " "imposing spatial constraints with a mask consisting of " "disconnected components will leading to clustering instability " @@ -384,18 +383,18 @@ def get_parser(): metavar="Cluster mask", default=None, nargs="+", - help="(Hyperparameter): Specify the path to a Nifti1Image mask file" + help="(metaparameter): Specify the path to a Nifti1Image mask file" " to constrained functional clustering. If specifying a list of " "paths to multiple cluster masks, separate them by space.\n", ) - # dMRI hyperparameters + # dMRI metaparameters parser.add_argument( "-ml", metavar="Minimum fiber length for tracking", default=20, nargs="+", - help="(Hyperparameter): Include this flag to manually specify a " + help="(metaparameter): Include this flag to manually specify a " "minimum tract length (mm) for dmri connectome tracking. Default " "is 20. If you wish to iterate the pipeline across multiple " "minimums, separate the list by space (e.g. 10 30 50).\n", @@ -405,7 +404,7 @@ def get_parser(): metavar="Error margin", default=8, nargs="+", - help="(Hyperparameter): Distance (in the units of the streamlines, " + help="(metaparameter): Distance (in the units of the streamlines, " "usually mm). If any coordinate in the streamline is within this " "distance from the center of any voxel in the ROI, the filtering " "criterion is set to True for this streamline, otherwise False. " @@ -418,7 +417,7 @@ def get_parser(): default="det", nargs="+", choices=["det", "prob", "clos"], - help="(Hyperparameter): Include this flag to manually specify the " + help="(metaparameter): Include this flag to manually specify the " "statistical approach to tracking for dmri connectome estimation." " Options are: det (deterministic), closest (clos), and " "prob (probabilistic). Default is det. If you wish to iterate the" @@ -650,10 +649,14 @@ def build_workflow(args, retval): ) if environ.get('FSLDIR') is None: - raise EnvironmentError('FSLDIR not found! Be sure that this ' - 'environment variable is set and/or that FSL ' - 'has been properly installed before ' - 'proceeding.') + try: + raise EnvironmentError('FSLDIR not found! Be sure that this ' + 'environment variable is set and/or that ' + 'FSL has been properly installed before ' + 'proceeding.') + except EnvironmentError: + retval["return_code"] = 1 + return retval else: fsl_version = check_output('flirt -version | cut -f3 -d\" \"', shell=True).strip() @@ -661,8 +664,12 @@ def build_workflow(args, retval): print(f"{Fore.MAGENTA}FSL {fsl_version.decode()} detected: " f"FSLDIR={os.environ['FSLDIR']}") else: - raise EnvironmentError('Is your FSL installation corrupted? ' - 'Check permissions.') + try: + raise EnvironmentError('Is your FSL installation corrupted? ' + 'Check permissions.') + except EnvironmentError: + retval["return_code"] = 1 + return retval # Start timer now = datetime.datetime.now() timestamp = str(now.strftime("%Y-%m-%d %H:%M:%S")) @@ -744,9 +751,13 @@ def build_workflow(args, retval): and multi_subject_graph is None and multi_subject_multigraph is None ): - raise ValueError( - "\nError: You must include a subject ID in your command line call." - ) + try: + raise ValueError( + "\nYou must include a subject ID in your command line " + "call." + ) + except ValueError: + sys.exit(1) if ( multi_subject_graph or multi_subject_multigraph or graph or multi_graph @@ -754,7 +765,7 @@ def build_workflow(args, retval): if multi_subject_graph: if len(ID) != len(multi_subject_graph): print( - "Error: Length of ID list does not correspond to length of" + "\nLength of ID list does not correspond to length of" " input graph file list." ) retval["return_code"] = 1 @@ -762,7 +773,7 @@ def build_workflow(args, retval): if multi_subject_multigraph: if len(ID) != len(multi_subject_multigraph): print( - "Error: Length of ID list does not correspond to length of" + "\nLength of ID list does not correspond to length of" " input graph file list." ) retval["return_code"] = 1 @@ -770,7 +781,7 @@ def build_workflow(args, retval): if len(ID) > 1 and not multi_subject_graph and not \ multi_subject_multigraph: print( - "Error: Length of ID list does not correspond to length of" + "\nLength of ID list does not correspond to length of" " input graph file list." ) retval["return_code"] = 1 @@ -1168,9 +1179,13 @@ def build_workflow(args, retval): directget = None if (ID is None) and (func_file_list is None): - raise ValueError( - "\nError: You must include a subject ID in your command line call." - ) + try: + raise ValueError( + "\nYou must include a subject ID in your command line " + "call." + ) + except ValueError: + sys.exit(1) if func_file_list and isinstance(ID, list): if len(ID) != len(func_file_list): @@ -1855,19 +1870,26 @@ def build_workflow(args, retval): print(f"{Fore.BLUE}{', '.join(min_length_list)}") if error_margin: if float(roi_neighborhood_tol) <= float(error_margin): - raise ValueError( - 'roi_neighborhood_tol preset cannot be less than ' - 'the value of the structural connectome error_margin' - ' parameter.') + try: + raise ValueError( + 'roi_neighborhood_tol preset cannot be less than ' + 'the value of the structural connectome error_margin' + ' parameter.') + except ValueError: + sys.exit(1) + print(f"{Fore.GREEN}Using {Fore.BLUE}{error_margin}" f"mm{Fore.GREEN} error margin...") else: for em in error_margin_list: if float(roi_neighborhood_tol) <= float(em): - raise ValueError( - 'roi_neighborhood_tol preset cannot be less than ' - 'the value of the structural connectome error_margin' - ' parameter.') + try: + raise ValueError( + 'roi_neighborhood_tol preset cannot be less than ' + 'the value of the structural connectome ' + 'error_margin parameter.') + except ValueError: + sys.exit(1) print(f"{Fore.GREEN}Iterating ROI-streamline intersection " f"tolerance:") print(f"{Fore.BLUE}{', '.join(error_margin_list)}") @@ -1891,44 +1913,67 @@ def build_workflow(args, retval): print(f"{Fore.GREEN}Diffusion-Weighted Image:{Fore.BLUE}\n " f"{_dwi_file}") if not os.path.isfile(_dwi_file): - raise FileNotFoundError(f"{_dwi_file} does not exist. " - f"Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{_dwi_file} does not exist. " + f"Ensure " + f"that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) + print(f"{Fore.GREEN}B-Values:\n{Fore.BLUE} {_fbval}") print(f"{Fore.GREEN}B-Vectors:\n{Fore.BLUE} {_fbvec}") if not os.path.isfile(fbvec): - raise FileNotFoundError(f"{_fbvec} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{_fbvec} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) + if not os.path.isfile(fbval): - raise FileNotFoundError(f"{_fbval} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{_fbval} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) else: print(f"{Fore.GREEN}Diffusion-Weighted Image:\n " f"{Fore.BLUE}{dwi_file}") if not os.path.isfile(dwi_file): - raise FileNotFoundError(f"{dwi_file} does not exist. " - f"Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{dwi_file} does not exist. " + f"Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) print(f"{Fore.GREEN}B-Values:\n {Fore.BLUE}{fbval}") print(f"{Fore.GREEN}B-Vectors:\n {Fore.BLUE}{fbvec}") if not os.path.isfile(fbvec): - raise FileNotFoundError(f"{fbvec} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{fbvec} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) if not os.path.isfile(fbval): - raise FileNotFoundError(f"{fbval} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{fbval} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) if waymask is not None: print(f"{Fore.GREEN}Waymask:\n {Fore.BLUE}{waymask}") if not os.path.isfile(waymask): - raise FileNotFoundError(f"{waymask} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{waymask} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) conf = None k = None clust_mask = None @@ -1949,31 +1994,42 @@ def build_workflow(args, retval): for _func_file in func_file_list: print(f"{Fore.GREEN}BOLD Image:\n {Fore.BLUE}{_func_file}") if not os.path.isfile(_func_file): - raise FileNotFoundError( - f"{_func_file} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError( + f"{_func_file} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) else: print(f"{Fore.GREEN}BOLD Image:\n {Fore.BLUE}{func_file}") if not os.path.isfile(func_file): - raise FileNotFoundError(f"{func_file} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") - + try: + raise FileNotFoundError(f"{func_file} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) if conf_list: for _conf in conf_list: print(f"{Fore.GREEN}BOLD Confound Regressors:\n " f"{Fore.BLUE}{_conf}") if not os.path.isfile(_conf): - raise FileNotFoundError(f"{_conf} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{_conf} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) elif conf: print(f"{Fore.GREEN}BOLD Confound Regressors:\n {Fore.BLUE}{conf}") if not os.path.isfile(conf): - raise FileNotFoundError(f"{conf} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{conf} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) multimodal = False elif (func_file or func_file_list) and (dwi_file or dwi_file_list): multimodal = True @@ -1983,30 +2039,42 @@ def build_workflow(args, retval): for _func_file in func_file_list: print(f"{Fore.GREEN}BOLD Image:\n {Fore.BLUE}{_func_file}") if not os.path.isfile(_func_file): - raise FileNotFoundError( - f"{_func_file} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError( + f"{_func_file} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) else: print(f"{Fore.GREEN}BOLD Image:\n {Fore.BLUE}{func_file}") if not os.path.isfile(func_file): - raise FileNotFoundError(f"{func_file} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{func_file} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) if conf_list: for _conf in conf_list: print(f"{Fore.GREEN}BOLD Confound Regressors:\n " f"{Fore.BLUE}{_conf}") if not os.path.isfile(_conf): - raise FileNotFoundError(f"{_conf} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{_conf} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) elif conf: print(f"{Fore.GREEN}BOLD Confound Regressors:\n {Fore.BLUE}{conf}") if not os.path.isfile(conf): - raise FileNotFoundError(f"{conf} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{conf} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) if dwi_file_list: for (_dwi_file, _fbval, _fbvec, _anat_file) in list( zip(dwi_file_list, fbval_list, fbvec_list, anat_file_list) @@ -2014,85 +2082,120 @@ def build_workflow(args, retval): print(f"{Fore.GREEN}Diffusion-Weighted Image:\n " f"{Fore.BLUE}{_dwi_file}") if not os.path.isfile(_dwi_file): - raise FileNotFoundError(f"{_dwi_file} does not exist. " - f"Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{_dwi_file} does not exist." + f" Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) print(f"{Fore.GREEN}B-Values:\n {Fore.BLUE}{_fbval}") print(f"{Fore.GREEN}B-Vectors:\n {Fore.BLUE}{_fbvec}") if not os.path.isfile(_fbvec): - raise FileNotFoundError(f"{_fbvec} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{_fbvec} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) if not os.path.isfile(_fbval): - raise FileNotFoundError(f"{_fbval} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{_fbval} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) else: print(f"{Fore.GREEN}Diffusion-Weighted Image:\n " f"{Fore.BLUE}{dwi_file}") if not os.path.isfile(dwi_file): - raise FileNotFoundError(f"{dwi_file} does not exist. " - f"Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{dwi_file} does not exist. " + f"Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) print(f"{Fore.GREEN}B-Values:\n {Fore.BLUE}{fbval}") print(f"{Fore.GREEN}B-Vectors:\n {Fore.BLUE}{fbvec}") if not os.path.isfile(fbvec): - raise FileNotFoundError(f"{fbvec} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{fbvec} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) if not os.path.isfile(fbval): - raise FileNotFoundError(f"{fbval} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{fbval} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) if waymask is not None: print(f"{Fore.GREEN}Waymask:\n {Fore.BLUE}{waymask}") if not os.path.isfile(waymask): - raise FileNotFoundError(f"{waymask} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{waymask} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) else: multimodal = False if roi is not None and roi is not 'None': print(f"{Fore.GREEN}ROI:\n {Fore.BLUE}{roi}") if not os.path.isfile(roi): - raise FileNotFoundError(f"{roi} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{roi} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) if anat_file or anat_file_list: if anat_file_list and len(anat_file_list) > 1: for anat_file in anat_file_list: print(f"{Fore.GREEN}T1-Weighted Image:\n " f"{Fore.BLUE}{anat_file}") if not os.path.isfile(anat_file): - raise FileNotFoundError( - f"{anat_file} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError( + f"{anat_file} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) else: print(f"{Fore.GREEN}T1-Weighted Image:\n {Fore.BLUE}{anat_file}") if not os.path.isfile(anat_file): - raise FileNotFoundError(f"{anat_file} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{anat_file} does not exist. " + f"Ensure that you are only " + f"specifying absolute paths.") + except FileNotFoundError: + sys.exit(1) if mask or mask_list: if mask_list and len(mask_list) > 1: for mask in mask_list: print(f"{Fore.GREEN}Brain Mask Image:\n {Fore.BLUE}{mask}") if not os.path.isfile(mask): - raise FileNotFoundError( - f"{mask} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError( + f"{mask} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) else: print(f"{Fore.GREEN}Brain Mask Image:\n {Fore.BLUE}{mask}") if not os.path.isfile(mask): - raise FileNotFoundError(f"{mask} does not exist. Ensure " - f"that you are only specifying " - f"absolute paths.") + try: + raise FileNotFoundError(f"{mask} does not exist. Ensure " + f"that you are only specifying " + f"absolute paths.") + except FileNotFoundError: + sys.exit(1) print(Style.RESET_ALL) print( "\n-------------------------------------------------------------------" diff --git a/pynets/core/interfaces.py b/pynets/core/interfaces.py index b48129fd..1a46d8f4 100755 --- a/pynets/core/interfaces.py +++ b/pynets/core/interfaces.py @@ -157,7 +157,8 @@ def _run_interface(self, runtime): uatlas = f"{runtime.cwd}{self.inputs.atlas}{'.nii.gz'}" if self.inputs.clustering is False: [uatlas, - labels] = nodemaker.enforce_hem_distinct_consecutive_labels( + labels] = \ + nodemaker.enforce_hem_distinct_consecutive_labels( uatlas, label_names=labels) if self.inputs.parc is True: parcel_list = nodemaker.gen_img_list(uatlas) @@ -343,7 +344,8 @@ def _run_interface(self, runtime): print(f"Coordinates:\n{coords}") - print(f"Labels:\n{textwrap.shorten(str(labels), width=1000, placeholder='...')}") + print(f"Labels:\n" + f"{textwrap.shorten(str(labels), width=1000, placeholder='...')}") assert len(coords) == len(labels) @@ -652,8 +654,8 @@ def create_bs_imgs(ts_data, block_size, clust_mask_corr_img): import nibabel as nib from nilearn.masking import unmask from pynets.fmri.estimation import timeseries_bootstrap - boot_series = timeseries_bootstrap(ts_data.func, - block_size)[0].astype('float32') + boot_series = timeseries_bootstrap( + ts_data.func, block_size)[0].astype('float32') return unmask(boot_series, clust_mask_corr_img) def run_bs_iteration(i, ts_data, work_dir, local_corr, diff --git a/pynets/core/nodemaker.py b/pynets/core/nodemaker.py index 08ae646c..de1f34f7 100644 --- a/pynets/core/nodemaker.py +++ b/pynets/core/nodemaker.py @@ -393,7 +393,7 @@ def get_node_membership( parcel_list_res = gen_img_list(par_tmp) else: parcel_list_res = None - + # Determine whether input is from 17-networks or 7-networks seven_nets = [ "Vis", @@ -973,7 +973,8 @@ def gen_img_list(uatlas): par_max = len(bna_data_for_coords_uniq) - 1 img_stack = [] for idx in range(1, par_max + 1): - roi_img = bna_data.astype("uint16") == bna_data_for_coords_uniq[idx].astype("uint16") + roi_img = bna_data.astype("uint16") == \ + bna_data_for_coords_uniq[idx].astype("uint16") img_stack.append(roi_img.astype("uint16")) img_stack = np.array(img_stack) @@ -1658,8 +1659,8 @@ def node_gen(coords, parcel_list, labels, dir_path, ID, parc, atlas, uatlas): else: label_intensities = labels - [net_parcels_map_nifti, _] = nodemaker.create_parcel_atlas(parcel_list, - label_intensities) + [net_parcels_map_nifti, _] = \ + nodemaker.create_parcel_atlas(parcel_list, label_intensities) coords = list(tuple(x) for x in coords) diff --git a/pynets/core/thresholding.py b/pynets/core/thresholding.py index 794b42eb..718d6204 100644 --- a/pynets/core/thresholding.py +++ b/pynets/core/thresholding.py @@ -36,8 +36,8 @@ def threshold_absolute(W, thr, copy=True): References ---------- - .. [1] Complex network measures of brain connectivity: Uses and interpretations. - Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. + .. [1] Complex network measures of brain connectivity: Uses and + interpretations. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. """ if copy: @@ -134,8 +134,8 @@ def normalize(W): References ---------- - .. [1] Complex network measures of brain connectivity: Uses and interpretations. - Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. + .. [1] Complex network measures of brain connectivity: Uses and + interpretations. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. """ W /= np.max(np.abs(W)) @@ -158,8 +158,8 @@ def standardize(W): References ---------- - .. [1] Complex network measures of brain connectivity: Uses and interpretations. - Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. + .. [1] Complex network measures of brain connectivity: Uses and + interpretations. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. """ W = (W - np.min(W)) / np.ptp(W) @@ -177,9 +177,11 @@ def density_thresholding(conn_matrix, thr, max_iters=10000, interval=0.01): thr : float Density value between 0-1. max_iters : int - Maximum number of iterations for performing absolute thresholding. Default is 1000. + Maximum number of iterations for performing absolute thresholding. + Default is 1000. interval : float - Interval for increasing the absolute threshold for each iteration. Default is 0.01. + Interval for increasing the absolute threshold for each iteration. + Default is 0.01. Returns ------- @@ -192,8 +194,8 @@ def density_thresholding(conn_matrix, thr, max_iters=10000, interval=0.01): Comparing brain networks of different size and connectivity density using graph theory. PLoS ONE. https://doi.org/10.1371/journal.pone.0013701 - .. [2] Complex network measures of brain connectivity: Uses and interpretations. - Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. + .. [2] Complex network measures of brain connectivity: Uses and + interpretations. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. """ from pynets.core import thresholding @@ -230,7 +232,8 @@ def density_thresholding(conn_matrix, thr, max_iters=10000, interval=0.01): i = i + 1 else: print( - "Density of raw matrix is already greater than or equal to the target density requested" + "Density of raw matrix is already greater than or equal to the " + "target density requested" ) return conn_matrix @@ -271,8 +274,8 @@ def thr2prob(W, copy=True): References ---------- - .. [1] Complex network measures of brain connectivity: Uses and interpretations. - Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. + .. [1] Complex network measures of brain connectivity: Uses and + interpretations. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. """ if copy: @@ -301,8 +304,8 @@ def binarize(W, copy=True): References ---------- - .. [1] Complex network measures of brain connectivity: Uses and interpretations. - Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. + .. [1] Complex network measures of brain connectivity: Uses and + interpretations. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. """ if copy: @@ -333,8 +336,8 @@ def invert(W, copy=False): References ---------- - .. [1] Complex network measures of brain connectivity: Uses and interpretations. - Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. + .. [1] Complex network measures of brain connectivity: Uses and + interpretations. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. """ if copy: @@ -391,8 +394,8 @@ def weight_conversion(W, wcm, copy=True): References ---------- - .. [1] Complex network measures of brain connectivity: Uses and interpretations. - Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. + .. [1] Complex network measures of brain connectivity: Uses and + interpretations. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. """ if wcm == "binarize": @@ -421,8 +424,8 @@ def autofix(W, copy=True): References ---------- - .. [1] Complex network measures of brain connectivity: Uses and interpretations. - Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. + .. [1] Complex network measures of brain connectivity: Uses and + interpretations. Rubinov M, Sporns O (2010) NeuroImage 52:1059-69. """ if copy: @@ -446,7 +449,8 @@ def autofix(W, copy=True): def disparity_filter(G, weight="weight"): """ - Compute significance scores (alpha) for weighted edges in G as defined in Serrano et al. 2009. + Compute significance scores (alpha) for weighted edges in G as defined in + Serrano et al. 2009. Parameters ---------- @@ -459,12 +463,13 @@ def disparity_filter(G, weight="weight"): Returns ------- B : Object - Weighted NetworkX graph with a significance score (alpha) assigned to each edge. + Weighted NetworkX graph with a significance score (alpha) assigned to + each edge. References ---------- - .. [1] M. A. Serrano et al. (2009) Extracting the Multiscale backbone of complex weighted networks. - PNAS, 106:16, pp. 6483-6488. + .. [1] M. A. Serrano et al. (2009) Extracting the Multiscale backbone of + complex weighted networks. PNAS, 106:16, pp. 6483-6488. """ from scipy import integrate @@ -485,9 +490,8 @@ def disparity_filter(G, weight="weight"): alpha_ij_out = ( 1 - (k_out - 1) - * integrate.quad(lambda x: (1 - x) ** (k_out - 2), 0, p_ij_out)[ - 0 - ] + * integrate.quad(lambda x: (1 - x) ** (k_out - 2), 0, + p_ij_out)[0] ) N.add_edge( u, v, weight=w, alpha_out=float(f"{alpha_ij_out:.4f}")) @@ -510,7 +514,8 @@ def disparity_filter(G, weight="weight"): alpha_ij_in = ( 1 - (k_in - 1) - * integrate.quad(lambda x: (1 - x) ** (k_in - 2), 0, p_ij_in)[0] + * integrate.quad(lambda x: (1 - x) ** (k_in - 2), 0, + p_ij_in)[0] ) N.add_edge( v, u, weight=w, alpha_in=float(f"{alpha_ij_in:.4f}")) @@ -528,7 +533,8 @@ def disparity_filter(G, weight="weight"): alpha_ij = ( 1 - (k - 1) - * integrate.quad(lambda x: (1 - x) ** (k - 2), 0, p_ij)[0] + * integrate.quad(lambda x: (1 - x) ** (k - 2), 0, + p_ij)[0] ) B.add_edge(u, v, weight=w, alpha=float(f"{alpha_ij:.4f}")) else: @@ -538,7 +544,8 @@ def disparity_filter(G, weight="weight"): def disparity_filter_alpha_cut(G, weight="weight", alpha_t=0.4, cut_mode="or"): """ - Compute significance scores (alpha) for weighted edges in G as defined in Serrano et al. 2009. + Compute significance scores (alpha) for weighted edges in G as defined in + Serrano et al. 2009. Parameters ---------- @@ -553,32 +560,36 @@ def disparity_filter_alpha_cut(G, weight="weight", alpha_t=0.4, cut_mode="or"): used to select the surviving edges. Default is 0.4. cut_mode : str - In the case of directed graphs. It represents the logic operation to filter out edges - that do not pass the threshold value, combining the alpha_in and alpha_out attributes - resulting from the disparity_filter function. Default is 'or'. Possible strings: 'or', 'and'. + In the case of directed graphs. It represents the logic operation + to filter out edges that do not pass the threshold value, + combining the alpha_in and alpha_out attributes resulting from + the disparity_filter function. Default is 'or'. + Possible strings: 'or', 'and'. + Returns ------- B : Object - Weighted NetworkX graph with a significance score (alpha) assigned to each edge. The resulting - graph contains only edges that survived from the filtering with the alpha_t threshold. + Weighted NetworkX graph with a significance score (alpha) assigned to + each edge. The resulting graph contains only edges that survived from + the filtering with the alpha_t threshold. References ---------- - .. [1] M. A. Serrano et al. (2009) Extracting the Multiscale backbone of complex weighted networks. - PNAS, 106:16, pp. 6483-6488. + .. [1] M. A. Serrano et al. (2009) Extracting the Multiscale backbone of + complex weighted networks. PNAS, 106:16, pp. 6483-6488. """ - if nx.is_directed(G): # Directed case: + if nx.is_directed(G): B = nx.DiGraph() for u, v, w in G.edges(data=True): try: alpha_in = w["alpha_in"] - except KeyError: # there is no alpha_in, so we assign 1. It will never pass the cut + except KeyError: alpha_in = 1 try: alpha_out = w["alpha_out"] - except KeyError: # there is no alpha_out, so we assign 1. It will never pass the cut + except KeyError: alpha_out = 1 if cut_mode == "or": @@ -594,7 +605,7 @@ def disparity_filter_alpha_cut(G, weight="weight", alpha_t=0.4, cut_mode="or"): for u, v, w in G.edges(data=True): try: alpha = w["alpha"] - except KeyError: # there is no alpha, so we assign 1. It will never pass the cut + except KeyError: alpha = 1 if alpha < alpha_t: @@ -606,10 +617,11 @@ def disparity_filter_alpha_cut(G, weight="weight", alpha_t=0.4, cut_mode="or"): def weight_to_distance(G): """ Inverts all the edge weights so they become equivalent to distance measure. - With a weight, the higher the value the stronger the connection. With a distance, - the higher the value the "weaker" the connection. In this case there is no - measurement unit for the distance, as it is just a conversion from the weights. - The distances can be accessed in each node's property with constants. + With a weight, the higher the value the stronger the connection. With a + distance, the higher the value the "weaker" the connection. In this case + there is no measurement unit for the distance, as it is just a conversion + from the weights. The distances can be accessed in each node's property + with constants. Parameters ---------- @@ -669,15 +681,17 @@ def knn(conn_matrix, k): def local_thresholding_prop(conn_matrix, thr): """ - Threshold the adjacency matrix by building from the minimum spanning tree (MST) and adding - successive N-nearest neighbour degree graphs to achieve target proportional threshold. + Threshold the adjacency matrix by building from the minimum spanning tree + (MST) and adding successive N-nearest neighbour degree graphs to achieve + target proportional threshold. Parameters ---------- conn_matrix : array Weighted NxN matrix. thr : float - A proportional threshold, between 0 and 1, to achieve through local thresholding. + A proportional threshold, between 0 and 1, to achieve through local + thresholding. Returns ------- @@ -686,10 +700,11 @@ def local_thresholding_prop(conn_matrix, thr): References ---------- - .. [1] Alexander-Bloch, A. F., Gogtay, N., Meunier, D., Birn, R., Clasen, L., - Lalonde, F., … Bullmore, E. T. (2010). Disrupted modularity and local - connectivity of brain functional networks in childhood-onset schizophrenia. - Frontiers in Systems Neuroscience. https://doi.org/10.3389/fnsys.2010.00147 + .. [1] Alexander-Bloch, A. F., Gogtay, N., Meunier, D., Birn, R., Clasen, + L., Lalonde, F., … Bullmore, E. T. (2010). Disrupted modularity and local + connectivity of brain functional networks in childhood-onset + schizophrenia. Frontiers in Systems Neuroscience. + https://doi.org/10.3389/fnsys.2010.00147 .. [2] Tewarie, P., van Dellen, E., Hillebrand, A., & Stam, C. J. (2015). The minimum spanning tree: An unbiased method for brain network analysis. NeuroImage. https://doi.org/10.1016/j.neuroimage.2014.10.015 @@ -717,7 +732,8 @@ def local_thresholding_prop(conn_matrix, thr): if len_edges > edgenum: print( - f"Warning: The minimum spanning tree already has: {len_edges} edges, select more edges. Local Threshold " + f"Warning: The minimum spanning tree already has: {len_edges} " + f"edges, select more edges. Local Threshold " f"will be applied by just retaining the Minimum Spanning Tree") conn_matrix_thr = nx.to_numpy_array(G) return conn_matrix_thr @@ -727,7 +743,8 @@ def local_thresholding_prop(conn_matrix, thr): while ( len_edges < edgenum and k <= np.shape(conn_matrix)[0] - and (len(len_edge_list[-fail_tol:]) - len(set(len_edge_list[-fail_tol:]))) + and (len(len_edge_list[-fail_tol:]) - + len(set(len_edge_list[-fail_tol:]))) < (fail_tol - 1) ): # print(k) @@ -772,7 +789,8 @@ def local_thresholding_prop(conn_matrix, thr): except ValueError: import sys - print("MST thresholding failed. Check raw graph output manually for debugging.") + print("MST thresholding failed. Check raw graph output manually for" + " debugging.") sys.exit(1) @@ -817,13 +835,11 @@ def perform_thresholding( edge_threshold = f"{str(thr_perc)}%" G1 = thresholding.disparity_filter( nx.from_numpy_array(np.abs(conn_matrix))) - # G2 = nx.Graph([(u, v, d) for u, v, d in G1.edges(data=True) if d['alpha'] < thr]) print(f"Computing edge disparity significance with alpha = {thr}") print( - f"Filtered graph: nodes = {G1.number_of_nodes()}, edges = {G1.number_of_edges()}" + f"Filtered graph: nodes = {G1.number_of_nodes()}, " + f"edges = {G1.number_of_edges()}" ) - # print(f'Backbone graph: nodes = {G2.number_of_nodes()}, edges = {G2.number_of_edges()}") - # print(G2.edges(data=True)) conn_matrix_bin = thresholding.binarize(nx.to_numpy_array( G1, nodelist=sorted(G1.nodes()), dtype=np.float64)) # Enforce original dimensionality by padding with zeros. @@ -884,23 +900,26 @@ def thresh_func( check_consistency=True, ): """ - Threshold a functional connectivity matrix using any of a variety of methods. + Threshold a functional connectivity matrix using any of a variety of + methods. Parameters ---------- dens_thresh : bool - Indicates whether a target graph density is to be used as the basis for thresholding. + Indicates whether a target graph density is to be used as the basis + for thresholding. thr : float - A value, between 0 and 1, to threshold the graph using any variety of methods - triggered through other options. + A value, between 0 and 1, to threshold the graph using any variety of + methods triggered through other options. conn_matrix : array Adjacency matrix stored as an m x n array of nodes and edges. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. ID : str A subject id or other unique identifier. dir_path : str @@ -908,13 +927,14 @@ def thresh_func( roi : str File path to binarized/boolean region-of-interest Nifti1Image file. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's. min_span_tree : bool Indicates whether local thresholding from the Minimum Spanning Tree should be used. smooth : int - Smoothing width (mm fwhm) to apply to time-series when extracting signal from ROI's. + Smoothing width (mm fwhm) to apply to time-series when extracting + signal from ROI's. disp_filt : bool Indicates whether local thresholding using a disparity filter and 'backbone network' should be used. @@ -939,7 +959,8 @@ def thresh_func( hpass : float High-pass filter values (Hz) to apply to node-extracted time-series. extract_strategy : str - The name of a valid function used to reduce the time-series region extraction. + The name of a valid function used to reduce the time-series region + extraction. Returns ------- @@ -948,23 +969,26 @@ def thresh_func( edge_threshold : str The string percentage representation of thr. est_path : str - File path to the thresholded graph, conn_matrix_thr, saved as a numpy array in .npy format. + File path to the thresholded graph, conn_matrix_thr, saved as a numpy + array in .npy format. thr : float - The value, between 0 and 1, used to threshold the graph using any variety of methods - triggered through other options. + The value, between 0 and 1, used to threshold the graph using any + variety of methods triggered through other options. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. roi : str File path to binarized/boolean region-of-interest Nifti1Image file. smooth : int - Smoothing width (mm fwhm) to apply to time-series when extracting signal from ROI's. + Smoothing width (mm fwhm) to apply to time-series when extracting + signal from ROI's. prune : bool Indicates whether to prune final graph of disconnected nodes/isolates. ID : str @@ -988,7 +1012,8 @@ def thresh_func( hpass : float High-pass filter values (Hz) to apply to node-extracted time-series. extract_strategy : str - The name of a valid function used to reduce the time-series region extraction. + The name of a valid function used to reduce the time-series region + extraction. References ---------- @@ -1029,7 +1054,8 @@ def thresh_func( ), ) - [thr_type, edge_threshold, conn_matrix_thr] = thresholding.perform_thresholding( + [thr_type, edge_threshold, conn_matrix_thr] = \ + thresholding.perform_thresholding( conn_matrix, thr, min_span_tree, dens_thresh, disp_filt) if not nx.is_connected(nx.from_numpy_matrix(conn_matrix_thr)): @@ -1111,7 +1137,8 @@ def thresh_struct( check_consistency=True, ): """ - Threshold a structural connectivity matrix using any of a variety of methods. + Threshold a structural connectivity matrix using any of a variety of + methods. Parameters ---------- @@ -1119,16 +1146,17 @@ def thresh_struct( Indicates whether a target graph density is to be used as the basis for thresholding. thr : float - A value, between 0 and 1, to threshold the graph using any variety of methods - triggered through other options. + A value, between 0 and 1, to threshold the graph using any variety of + methods triggered through other options. conn_matrix : array Adjacency matrix stored as an m x n array of nodes and edges. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. ID : str A subject id or other unique identifier. dir_path : str @@ -1136,8 +1164,8 @@ def thresh_struct( roi : str File path to binarized/boolean region-of-interest Nifti1Image file. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's. min_span_tree : bool Indicates whether local thresholding from the Minimum Spanning Tree should be used. @@ -1171,7 +1199,8 @@ def thresh_struct( streams : str File path to save streamline array sequence in .trk format. directget : str - The statistical approach to tracking. Options are: det (deterministic), closest (clos), boot (bootstrapped), + The statistical approach to tracking. Options are: + det (deterministic), closest (clos), boot (bootstrapped), and prob (probabilistic). min_length : int Minimum fiber length threshold in mm to restrict tracking. @@ -1183,19 +1212,21 @@ def thresh_struct( edge_threshold : str The string percentage representation of thr. est_path : str - File path to the thresholded graph, conn_matrix_thr, saved as a numpy array in .npy format. + File path to the thresholded graph, conn_matrix_thr, saved as a numpy + array in .npy format. thr : float - The value, between 0 and 1, used to threshold the graph using any variety of methods - triggered through other options. + The value, between 0 and 1, used to threshold the graph using any + variety of methods triggered through other options. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. roi : str File path to binarized/boolean region-of-interest Nifti1Image file. prune : bool @@ -1227,7 +1258,8 @@ def thresh_struct( streams : str File path to save streamline array sequence in .trk format. directget : str - The statistical approach to tracking. Options are: det (deterministic), closest (clos), boot (bootstrapped), + The statistical approach to tracking. Options are: + det (deterministic), closest (clos), boot (bootstrapped), and prob (probabilistic). min_length : int Minimum fiber length threshold in mm to restrict tracking. diff --git a/pynets/core/utils.py b/pynets/core/utils.py index 8d77b296..c429e399 100644 --- a/pynets/core/utils.py +++ b/pynets/core/utils.py @@ -51,7 +51,8 @@ def prune_suffices(res): def do_dir_path(atlas, outdir): """ - Creates an atlas subdirectory from the base directory of the given subject's input file. + Creates an atlas subdirectory from the base directory of the given + subject's input file. Parameters ---------- @@ -78,7 +79,8 @@ def do_dir_path(atlas, outdir): os.makedirs(dir_path, exist_ok=True) elif atlas is None: try: - raise ValueError("Error: cannot create directory for a null atlas!") + raise ValueError("Error: cannot create directory for a null " + "atlas!") except ValueError: import sys sys.exit(1) @@ -134,30 +136,33 @@ def create_est_path_func( extract_strategy, ): """ - Name the thresholded functional connectivity matrix file based on relevant graph-generating parameters. + Name the thresholded functional connectivity matrix file based on + relevant graph-generating parameters. Parameters ---------- ID : str A subject id or other unique identifier. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. thr : float - A value, between 0 and 1, to threshold the graph using any variety of methods - triggered through other options. + A value, between 0 and 1, to threshold the graph using any variety of + methods triggered through other options. roi : str File path to binarized/boolean region-of-interest Nifti1Image file. dir_path : str Path to directory containing subject derivative data for given run. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's. smooth : int - Smoothing width (mm fwhm) to apply to time-series when extracting signal from ROI's. + Smoothing width (mm fwhm) to apply to time-series when extracting + signal from ROI's. thr_type : str Type of thresholding performed (e.g. prop, abs, dens, mst, disp) hpass : bool @@ -165,12 +170,14 @@ def create_est_path_func( parc : bool Indicates whether to use parcels instead of coordinates as ROI nodes. extract_strategy : str - The name of a valid function used to reduce the time-series region extraction. + The name of a valid function used to reduce the time-series region + extraction. Returns ------- est_path : str - File path to .npy file containing graph with all specified combinations of hyperparameter characteristics. + File path to .npy file containing graph with all specified + combinations of hyperparameter characteristics. """ import os @@ -260,28 +267,30 @@ def create_est_path_diff( error_margin, ): """ - Name the thresholded structural connectivity matrix file based on relevant graph-generating parameters. + Name the thresholded structural connectivity matrix file based on + relevant graph-generating parameters. Parameters ---------- ID : str A subject id or other unique identifier. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for + partial correlation). sps type is used by default. thr : float - A value, between 0 and 1, to threshold the graph using any variety of methods - triggered through other options. + A value, between 0 and 1, to threshold the graph using any variety of + methods triggered through other options. roi : str File path to binarized/boolean region-of-interest Nifti1Image file. dir_path : str Path to directory containing subject derivative data for given run. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's. target_samples : int Total number of streamline samples specified to generate streams. track_type : str @@ -291,7 +300,8 @@ def create_est_path_diff( parc : bool Indicates whether to use parcels instead of coordinates as ROI nodes. directget : str - The statistical approach to tracking. Options are: det (deterministic), closest (clos), boot (bootstrapped), + The statistical approach to tracking. Options are: + det (deterministic), closest (clos), boot (bootstrapped), and prob (probabilistic). min_length : int Minimum fiber length threshold in mm to restrict tracking. @@ -380,38 +390,43 @@ def create_raw_path_func( extract_strategy, ): """ - Name the raw functional connectivity matrix file based on relevant graph-generating parameters. + Name the raw functional connectivity matrix file based on relevant + graph-generating parameters. Parameters ---------- ID : str A subject id or other unique identifier. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for + partial correlation). sps type is used by default. roi : str File path to binarized/boolean region-of-interest Nifti1Image file. dir_path : str Path to directory containing subject derivative data for given run. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's. smooth : int - Smoothing width (mm fwhm) to apply to time-series when extracting signal from ROI's. + Smoothing width (mm fwhm) to apply to time-series when extracting + signal from ROI's. hpass : bool High-pass filter values (Hz) to apply to node-extracted time-series. parc : bool Indicates whether to use parcels instead of coordinates as ROI nodes. extract_strategy : str - The name of a valid function used to reduce the time-series region extraction. + The name of a valid function used to reduce the time-series region + extraction. Returns ------- est_path : str - File path to .npy file containing graph with all specified combinations of hyperparameter characteristics. + File path to .npy file containing graph with all specified + combinations of hyperparameter characteristics. """ import os @@ -617,7 +632,8 @@ def create_csv_path(dir_path, est_path): if not os.path.isdir(namer_dir): os.makedirs(namer_dir, exist_ok=True) - out_path = f"{namer_dir}/metrics_{est_path.split('/')[-1].split('.npy')[0]}.csv" + out_path = f"{namer_dir}/metrics_" \ + f"{est_path.split('/')[-1].split('.npy')[0]}.csv" return out_path @@ -828,55 +844,66 @@ def pass_meta_outs( Parameters ---------- conn_model_iterlist : list - List of connectivity estimation model parameters (e.g. corr for correlation, cov for covariance, - sps for precision covariance, partcorr for partial correlation). sps type is used by default. + List of connectivity estimation model parameters (e.g. corr for + correlation, cov for covariance, sps for precision covariance, + partcorr for partial correlation). sps type is used by default. est_path_iterlist : list - List of file paths to .npy file containing graph with thresholding applied. + List of file paths to .npy file containing graph with thresholding + applied. network_iterlist : list - List of resting-state networks based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the - study of brain subgraphs. + List of resting-state networks based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. thr_iterlist : list - List of values, between 0 and 1, to threshold the graph using any variety of methods - triggered through other options. + List of values, between 0 and 1, to threshold the graph using any + variety of methods triggered through other options. prune_iterlist : list - List of booleans indicating whether final graphs were pruned of disconnected nodes/isolates. + List of booleans indicating whether final graphs were pruned of + disconnected nodes/isolates. ID_iterlist : list List of repeated subject id strings. roi_iterlist : list - List of file paths to binarized/boolean region-of-interest Nifti1Image files. + List of file paths to binarized/boolean region-of-interest + Nifti1Image files. norm_iterlist : list Indicates method of normalizing resulting graph. binary_iterlist : list - List of booleans indicating whether resulting graph edges to form an unweighted graph were binarized. + List of booleans indicating whether resulting graph edges to form an + unweighted graph were binarized. Returns ------- conn_model_iterlist : list - List of connectivity estimation model parameters (e.g. corr for correlation, cov for covariance, - sps for precision covariance, partcorr for partial correlation). sps type is used by default. + List of connectivity estimation model parameters (e.g. corr for + correlation, cov for covariance, sps for precision covariance, + partcorr for partial correlation). sps type is used by default. est_path_iterlist : list - List of file paths to .npy file containing graph with thresholding applied. + List of file paths to .npy file containing graph with thresholding + applied. network_iterlist : list - List of resting-state networks based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the - study of brain subgraphs. + List of resting-state networks based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. thr_iterlist : list - List of values, between 0 and 1, to threshold the graph using any variety of methods - triggered through other options. + List of values, between 0 and 1, to threshold the graph using + any variety of methods triggered through other options. prune_iterlist : list - List of booleans indicating whether final graphs were pruned of disconnected nodes/isolates. + List of booleans indicating whether final graphs were pruned of + disconnected nodes/isolates. ID_iterlist : list List of repeated subject id strings. roi_iterlist : list - List of file paths to binarized/boolean region-of-interest Nifti1Image files. + List of file paths to binarized/boolean region-of-interest + Nifti1Image files. norm_iterlist : list Indicates method of normalizing resulting graph. binary_iterlist : list - List of booleans indicating whether resulting graph edges to form an unweighted graph were binarized. + List of booleans indicating whether resulting graph edges to form an + unweighted graph were binarized. embed_iterlist : list - List of booleans indicating whether omnibus embedding of graph population was performed. + List of booleans indicating whether omnibus embedding of graph + population was performed. multimodal_iterlist : list - List of booleans indicating whether multiple modalities of input data have been specified. - + List of booleans indicating whether multiple modalities of input data + have been specified. """ return ( @@ -908,16 +935,17 @@ def pass_meta_ins( Parameters ---------- conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. est_path : str File path to .npy file containing graph with thresholding applied. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') + used to filter nodes in the study of brain subgraphs. thr : float - A value, between 0 and 1, to threshold the graph using any variety of methods - triggered through other options. + A value, between 0 and 1, to threshold the graph using any variety of + methods triggered through other options. prune : bool Indicates whether to prune final graph of disconnected nodes/isolates. ID : str @@ -933,16 +961,17 @@ def pass_meta_ins( Returns ------- conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. est_path : str File path to .npy file containing graph with thresholding applied. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') + used to filter nodes in the study of brain subgraphs. thr : float - A value, between 0 and 1, to threshold the graph using any variety of methods - triggered through other options. + A value, between 0 and 1, to threshold the graph using any variety of + methods triggered through other options. prune : bool Indicates whether to prune final graph of disconnected nodes/isolates. ID : str @@ -954,7 +983,6 @@ def pass_meta_ins( binary : bool Indicates whether to binarize resulting graph edges to form an unweighted graph. - """ est_path_iterlist = est_path conn_model_iterlist = conn_model @@ -1015,77 +1043,94 @@ def pass_meta_ins_multi( Parameters ---------- conn_model_func : str - Functional connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision - covariance, partcorr for partial correlation). sps type is used by default. + Functional connectivity estimation model (e.g. corr for correlation, cov + for covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. est_path_func : str - File path to .npy file containing functional graph with thresholding applied. + File path to .npy file containing functional graph with thresholding + applied. network_func : str - Functional resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the - study of brain subgraphs. + Functional resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. thr_func : float - A value, between 0 and 1, to threshold the functional graph using any variety of methods - triggered through other options. + A value, between 0 and 1, to threshold the functional graph using any + variety of methods triggered through other options. prune_func : bool - Indicates whether to prune final functional graph of disconnected nodes/isolates. + Indicates whether to prune final functional graph of disconnected + nodes/isolates. ID_func : str A subject id or other unique identifier for the functional workflow. roi_func : str - File path to binarized/boolean region-of-interest Nifti1Image file applied to the functional data. + File path to binarized/boolean region-of-interest Nifti1Image file + applied to the functional data. norm_func : int Indicates method of normalizing resulting functional graph. binary_func : bool - Indicates whether to binarize resulting graph edges to form an unweighted functional graph. + Indicates whether to binarize resulting graph edges to form an + unweighted functional graph. conn_model_struct : str - Diffusion structural connectivity estimation model (e.g. corr for correlation, cov for covariance, - sps for precision covariance, partcorr for partial correlation). sps type is used by default. + Diffusion structural connectivity estimation model (e.g. corr for + correlation, cov for covariance, sps for precision covariance, partcorr + for partial correlation). sps type is used by default. est_path_struct : str - File path to .npy file containing diffusion structural graph with thresholding applied. + File path to .npy file containing diffusion structural graph with + thresholding applied. network_struct : str - Diffusion structural resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter - nodes in the study of brain subgraphs. + Diffusion structural resting-state network based on Yeo-7 and Yeo-17 + naming (e.g. 'Default') used to filter nodes in the study of brain + subgraphs. thr_struct : float - A value, between 0 and 1, to threshold the diffusion structural graph using any variety of methods - triggered through other options. + A value, between 0 and 1, to threshold the diffusion structural graph + using any variety of methods triggered through other options. prune_struct : bool - Indicates whether to prune final diffusion structural graph of disconnected nodes/isolates. + Indicates whether to prune final diffusion structural graph of + disconnected nodes/isolates. ID_struct : str - A subject id or other unique identifier for the diffusion structural workflow. + A subject id or other unique identifier for the diffusion structural + workflow. roi_struct : str - File path to binarized/boolean region-of-interest Nifti1Image file applied too the dwi data. + File path to binarized/boolean region-of-interest Nifti1Image file + applied too the dwi data. norm_struct : int Indicates method of normalizing resulting diffusion structural graph. binary_struct : bool - Indicates whether to binarize resulting diffusion structural graph edges to form an unweighted graph. - + Indicates whether to binarize resulting diffusion structural graph + edges to form an unweighted graph. Returns ------- conn_model_iterlist : list - List of connectivity estimation model parameters (e.g. corr for correlation, cov for covariance, - sps for precision covariance, partcorr for partial correlation). sps type is used by default. + List of connectivity estimation model parameters (e.g. corr for + correlation, cov for covariance, sps for precision covariance, partcorr + for partial correlation). sps type is used by default. est_path_iterlist : list - List of file paths to .npy file containing graph with thresholding applied. + List of file paths to .npy file containing graph with thresholding + applied. network_iterlist : list - List of resting-state networks based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the - study of brain subgraphs. + List of resting-state networks based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. thr_iterlist : list - List of values, between 0 and 1, to threshold the graph using any variety of methods - triggered through other options. + List of values, between 0 and 1, to threshold the graph using any + variety of methods triggered through other options. prune_iterlist : list - List of booleans indicating whether final graphs were pruned of disconnected nodes/isolates. + List of booleans indicating whether final graphs were pruned of + disconnected nodes/isolates. ID_iterlist : list List of repeated subject id strings. roi_iterlist : list - List of file paths to binarized/boolean region-of-interest Nifti1Image files. + List of file paths to binarized/boolean region-of-interest + Nifti1Image files. norm_iterlist : list Indicates method of normalizing resulting graph. binary_iterlist : list - List of booleans indicating whether resulting graph edges to form an unweighted graph were binarized. + List of booleans indicating whether resulting graph edges to form an + unweighted graph were binarized. embed_iterlist : list - List of booleans indicating whether omnibus embedding of graph population was performed. + List of booleans indicating whether omnibus embedding of graph + population was performed. multimodal_iterlist : list - List of booleans indicating whether multiple modalities of input data have been specified. - + List of booleans indicating whether multiple modalities of input data + have been specified. """ est_path_iterlist = [est_path_func, est_path_struct] conn_model_iterlist = [conn_model_func, conn_model_struct] @@ -1352,13 +1397,17 @@ def check_est_path_existence(est_path_list): Parameters ---------- est_path_list : list - List of file paths to .npy file containing graph with thresholding applied. + List of file paths to .npy file containing graph with thresholding + applied. + Returns ------- est_path_list_ex : list - List of existing file paths to .npy file containing graph with thresholding applied. + List of existing file paths to .npy file containing graph with + thresholding applied. bad_ixs : int - List of indices in est_path_list with non-existent and/or corrupt files. + List of indices in est_path_list with non-existent and/or corrupt + files. """ est_path_list_ex = [] @@ -1390,8 +1439,8 @@ def save_coords_and_labels_to_pickle(coords, labels, dir_path, network): dir_path : str Path to directory containing subject derivative data for given run. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. Returns ------- @@ -1469,9 +1518,11 @@ def get_template_tf(template_name, vox_size): return template, template_mask, templateflow_home -def save_nifti_parcels_map(ID, dir_path, network, net_parcels_map_nifti, vox_size): +def save_nifti_parcels_map(ID, dir_path, network, net_parcels_map_nifti, + vox_size): """ - This function takes a Nifti1Image parcellation object resulting from some form of masking and saves it to disk. + This function takes a Nifti1Image parcellation object resulting from some + form of masking and saves it to disk. Parameters ---------- @@ -1480,19 +1531,19 @@ def save_nifti_parcels_map(ID, dir_path, network, net_parcels_map_nifti, vox_siz dir_path : str Path to directory containing subject derivative data for given run. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. net_parcels_map_nifti : Nifti1Image - A nibabel-based nifti image consisting of a 3D array with integer voxel intensities corresponding to ROI - membership. + A nibabel-based nifti image consisting of a 3D array with integer + voxel intensities corresponding to ROI membership. vox_size : str Voxel size in mm. (e.g. 2mm). Returns ------- net_parcels_nii_path : str - File path to Nifti1Image consisting of a 3D array with integer voxel intensities corresponding to ROI - membership. + File path to Nifti1Image consisting of a 3D array with integer voxel + intensities corresponding to ROI membership. """ import os @@ -1565,29 +1616,33 @@ def save_ts_to_file( roi : str File path to binarized/boolean region-of-interest Nifti1Image file. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. ID : str A subject id or other unique identifier. dir_path : str Path to directory containing subject derivative data for given run. ts_within_nodes : array - 2D m x n array consisting of the time-series signal for each ROI node where m = number of scans and - n = number of ROI's, where ROI's are parcel volumes. + 2D m x n array consisting of the time-series signal for each ROI node + where m = number of scans and n = number of ROI's, where ROI's are + parcel volumes. smooth : int - Smoothing width (mm fwhm) to apply to time-series when extracting signal from ROI's. + Smoothing width (mm fwhm) to apply to time-series when extracting + signal from ROI's. hpass : bool High-pass filter values (Hz) to apply to node-extracted time-series. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's for time-series extraction. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's for time-series extraction. extract_strategy : str - The name of a valid function used to reduce the time-series region extraction. + The name of a valid function used to reduce the time-series region + extraction. Returns ------- out_path_ts : str - Path to .npy file containing array of fMRI time-series extracted from nodes. + Path to .npy file containing array of fMRI time-series extracted from + nodes. """ import os diff --git a/pynets/dmri/track.py b/pynets/dmri/track.py index b8d7583c..d85202f0 100644 --- a/pynets/dmri/track.py +++ b/pynets/dmri/track.py @@ -431,12 +431,18 @@ def track_ensemble( ) as stream: hardcoded_params = yaml.load(stream) nthreads = hardcoded_params["nthreads"][0] - n_seeds_per_iter = hardcoded_params['tracking']["n_seeds_per_iter"][0] - max_length = hardcoded_params['tracking']["max_length"][0] - pft_back_tracking_dist = hardcoded_params['tracking']["pft_back_tracking_dist"][0] - pft_front_tracking_dist = hardcoded_params['tracking']["pft_front_tracking_dist"][0] - particle_count = hardcoded_params['tracking']["particle_count"][0] - min_separation_angle = hardcoded_params['tracking']["min_separation_angle"][0] + n_seeds_per_iter = \ + hardcoded_params['tracking']["n_seeds_per_iter"][0] + max_length = \ + hardcoded_params['tracking']["max_length"][0] + pft_back_tracking_dist = \ + hardcoded_params['tracking']["pft_back_tracking_dist"][0] + pft_front_tracking_dist = \ + hardcoded_params['tracking']["pft_front_tracking_dist"][0] + particle_count = \ + hardcoded_params['tracking']["particle_count"][0] + min_separation_angle = \ + hardcoded_params['tracking']["min_separation_angle"][0] stream.close() all_combs = list(itertools.product(step_list, curv_thr_list)) @@ -465,6 +471,8 @@ def track_ensemble( out_streams = [i for i in out_streams if i is not None and len(i) > 0] + out_streams = concatenate(out_streams, axis=0) + if len(out_streams) < 100: ix += 1 print("Fewer than 100 streamlines tracked on last iteration." @@ -474,10 +482,7 @@ def track_ensemble( roi_neighborhood_tol = float(roi_neighborhood_tol) * 1.05 min_length = float(min_length) * 0.95 else: - out_streams = concatenate(out_streams, axis=0) - - if float(ix) > len(all_combs): - break + ix -= 1 # Append streamline generators to prevent exponential growth # in memory consumption @@ -497,6 +502,9 @@ def track_ensemble( gc.collect() print(Style.RESET_ALL) + if float(ix) > len(all_combs): + break + if ix >= len(all_combs) and float(stream_counter) < float(target_samples): print(f"Tractography failed. >{len(all_combs)} consecutive sampling " f"iterations with <100 streamlines. Are you using a waymask? " @@ -803,8 +811,8 @@ def run_tracking(step_curv_combinations, atlas_data_wm_gm_int, recon_path, return None if waymask is not None and os.path.isfile(waymask_tmp_path): - waymask_data = np.asarray(nib.load(waymask_tmp_path).dataobj).astype( - "bool") + waymask_data = np.asarray(nib.load(waymask_tmp_path).dataobj + ).astype("bool") try: roi_proximal_streamlines = roi_proximal_streamlines[ utils.near_roi( @@ -812,7 +820,7 @@ def run_tracking(step_curv_combinations, atlas_data_wm_gm_int, recon_path, np.eye(4), waymask_data, tol=0, - mode="any", + mode="all", ) ] print("%s%s" % ("Waymask proximity: ", diff --git a/pynets/fmri/clustools.py b/pynets/fmri/clustools.py index 5bd4bdba..b43a23a1 100644 --- a/pynets/fmri/clustools.py +++ b/pynets/fmri/clustools.py @@ -73,18 +73,19 @@ def ncut(W, nbEigenValues): This function performs the first step of normalized cut spectral clustering. The normalized LaPlacian is calculated on the similarity matrix W, and top nbEigenValues eigenvectors are calculated. The number of eigenvectors - corresponds to the maximum number of classes (K) that will be produced by the - clustering algorithm. + corresponds to the maximum number of classes (K) that will be produced by + the clustering algorithm. Parameters ---------- W : array - Numpy array containing a symmetric #feature x #feature sparse matrix representing the - similarity between voxels, traditionally this matrix should be positive semidefinite, - but regularization is employed to allow negative matrix entries (Yu 2001). + Numpy array containing a symmetric #feature x #feature sparse matrix + representing the similarity between voxels, traditionally this matrix + should be positive semidefinite, but regularization is employed to + allow negative matrix entries (Yu 2001). nbEigenValues : int - Number of eigenvectors that should be calculated, this determines the maximum number - of clusters (K) that can be derived from the result. + Number of eigenvectors that should be calculated, this determines the + maximum number of clusters (K) that can be derived from the result. Returns ------- @@ -95,14 +96,16 @@ def ncut(W, nbEigenValues): References ---------- - .. [1] Stella Yu and Jianbo Shi, "Understanding Popout through Repulsion," Computer - Vision and Pattern Recognition, December, 2001. - .. [2] Shi, J., & Malik, J. (2000). Normalized cuts and image segmentation. IEEE - Transactions on Pattern Analysis and Machine Intelligence, 22(8), 888-905. - doi: 10.1109/34.868688. - .. [3] Yu, S. X., & Shi, J. (2003). Multiclass spectral clustering. Proceedings Ninth - IEEE International Conference on Computer Vision, (1), 313-319 vol.1. Ieee. - doi: 10.1109/ICCV.2003.1238361 + .. [1] Stella Yu and Jianbo Shi, "Understanding Popout through Repulsion," + Computer Vision and Pattern Recognition, December, 2001. + .. [2] Shi, J., & Malik, J. (2000). Normalized cuts and image + segmentation. IEEE + Transactions on Pattern Analysis and Machine Intelligence, 22(8), + 888-905. doi: 10.1109/34.868688. + .. [3] Yu, S. X., & Shi, J. (2003). Multiclass spectral clustering. + Proceedings Ninth + IEEE International Conference on Computer Vision, (1), 313-319 vol.1. + Ieee. doi: 10.1109/ICCV.2003.1238361 """ from scipy.sparse.linalg import eigsh @@ -151,21 +154,22 @@ def ncut(W, nbEigenValues): def discretisation(eigen_vec): """ This function performs the second step of normalized cut clustering which - assigns features to clusters based on the eigen vectors from the LaPlacian of - a similarity matrix. There are a few different ways to perform this task. Shi - and Malik (2000) iteratively bisect the features based on the positive and - negative loadings of the eigenvectors. Ng, Jordan and Weiss (2001) proposed to - perform K-means clustering on the rows of the eigenvectors. The method + assigns features to clusters based on the eigen vectors from the LaPlacian + of a similarity matrix. There are a few different ways to perform this + task. Shi and Malik (2000) iteratively bisect the features based on the + positive and negative loadings of the eigenvectors. Ng, Jordan and Weiss + (2001) proposed to perform K-means clustering on the rows of the + eigenvectors. The method implemented here was proposed by Yu and Shi (2003) and it finds a discrete solution by iteratively rotating a binarised set of vectors until they are - maximally similar to the the eigenvectors. An advantage of this method over K-means - is that it is _more_ deterministic, i.e. you should get very similar results - every time you run the algorithm on the same data. + maximally similar to the the eigenvectors. An advantage of this method + over K-means is that it is _more_ deterministic, i.e. you should get very + similar results every time you run the algorithm on the same data. - The number of clusters that the features are clustered into is determined by - the number of eignevectors (number of columns) in the input array eigen_vec. A - caveat of this method, is that number of resulting clusters is bound by the - number of eignevectors, but it may contain less. + The number of clusters that the features are clustered into is determined + by the number of eignevectors (number of columns) in the input array + eigen_vec. A caveat of this method, is that number of resulting clusters + is bound by the number of eignevectors, but it may contain less. Parameters ---------- @@ -184,14 +188,15 @@ def discretisation(eigen_vec): References ---------- - .. [1] Stella Yu and Jianbo Shi, "Understanding Popout through Repulsion," Computer - Vision and Pattern Recognition, December, 2001. - .. [2] Shi, J., & Malik, J. (2000). Normalized cuts and image segmentation. IEEE - Transactions on Pattern Analysis and Machine Intelligence, 22(8), 888-905. - doi: 10.1109/34.868688. - .. [3] Yu, S. X., & Shi, J. (2003). Multiclass spectral clustering. Proceedings Ninth - IEEE International Conference on Computer Vision, (1), 313-319 vol.1. Ieee. - doi: 10.1109/ICCV.2003.1238361 + .. [1] Stella Yu and Jianbo Shi, "Understanding Popout through Repulsion," + Computer Vision and Pattern Recognition, December, 2001. + .. [2] Shi, J., & Malik, J. (2000). Normalized cuts and image + segmentation. IEEE + Transactions on Pattern Analysis and Machine Intelligence, 22(8), + 888-905. doi: 10.1109/34.868688. + .. [3] Yu, S. X., & Shi, J. (2003). Multiclass spectral clustering. + Proceedings Ninth IEEE International Conference on Computer Vision, + (1), 313-319 vol.1. Ieee. doi: 10.1109/ICCV.2003.1238361 """ import scipy as sp @@ -232,8 +237,10 @@ def discretisation(eigen_vec): nbIterationsDiscretisation = 0 nbIterationsDiscretisationMax = 20 - # Iteratively rotate the discretised eigenvectors until they are maximally similar to the input eignevectors, - # this converges when the differences between the current solution and the previous solution differs by less + # Iteratively rotate the discretised eigenvectors until they are + # maximally similar to the input eignevectors, + # this converges when the differences between the current solution + # and the previous solution differs by less # than eps or we have reached the maximum number of itarations while exitLoop == 0: nbIterationsDiscretisation = nbIterationsDiscretisation + 1 @@ -293,12 +300,14 @@ def parcellate_ncut(W, k, mask_img): Parameters ---------- W : Compressed Sparse Matrix - A Scipy sparse matrix, with weights corresponding to the temporal/spatial - correlation between the time series from voxel i and voxel j. + A Scipy sparse matrix, with weights corresponding to the + temporal/spatial correlation between the time series from voxel i + and voxel j. k : int Numbers of clusters that will be generated. mask_img : Nifti1Image - 3D NIFTI file containing a mask, which restricts the voxels used in the analysis. + 3D NIFTI file containing a mask, which restricts the voxels used in + the analysis. References ---------- @@ -744,30 +753,34 @@ def __init__( k : int Numbers of clusters that will be generated. clust_type : str - Type of clustering to be performed (e.g. 'ward', 'kmeans', 'complete', 'average'). + Type of clustering to be performed (e.g. 'ward', 'kmeans', + 'complete', 'average'). local_corr : str - Type of local connectivity to use as the basis for clustering methods. Options are tcorr or scorr. - Default is tcorr. + Type of local connectivity to use as the basis for clustering + methods. Options are tcorr or scorr. Default is tcorr. outdir : str Path to base derivatives directory. conf : str - File path to a confound regressor file for reduce noise in the time-series when extracting from ROI's. + File path to a confound regressor file for reduce noise in the + time-series when extracting from ROI's. mask : str File path to a 3D NIFTI file containing a mask, which restricts the voxels used in the analysis. References ---------- - .. [1] Thirion, B., Varoquaux, G., Dohmatob, E., & Poline, J. B. (2014). - Which fMRI clustering gives good brain parcellations? + .. [1] Thirion, B., Varoquaux, G., Dohmatob, E., & Poline, J. B. + (2014). Which fMRI clustering gives good brain parcellations? Frontiers in Neuroscience. https://doi.org/10.3389/fnins.2014.00167 - .. [2] Bellec, P., Rosa-Neto, P., Lyttelton, O. C., Benali, H., & Evans, A. C. (2010). - Multi-level bootstrap analysis of stable clusters in resting-state fMRI. - NeuroImage. https://doi.org/10.1016/j.neuroimage.2010.02.082 - .. [3] Garcia-Garcia, M., Nikolaidis, A., Bellec, P., Craddock, R. C., Cheung, B., - Castellanos, F. X., & Milham, M. P. (2018). Detecting stable individual - differences in the functional organization of the human basal ganglia. - NeuroImage. https://doi.org/10.1016/j.neuroimage.2017.07.029 + .. [2] Bellec, P., Rosa-Neto, P., Lyttelton, O. C., Benali, H., & + Evans, A. C. (2010). Multi-level bootstrap analysis of stable + clusters in resting-state fMRI. NeuroImage. + https://doi.org/10.1016/j.neuroimage.2010.02.082 + .. [3] Garcia-Garcia, M., Nikolaidis, A., Bellec, P., + Craddock, R. C., Cheung, B., Castellanos, F. X., & Milham, M. P. + (2018). Detecting stable individual differences in the functional + organization of the human basal ganglia. NeuroImage. + https://doi.org/10.1016/j.neuroimage.2017.07.029 """ self.func_file = func_file diff --git a/pynets/fmri/estimation.py b/pynets/fmri/estimation.py index a50e68b3..10cc4cee 100644 --- a/pynets/fmri/estimation.py +++ b/pynets/fmri/estimation.py @@ -100,24 +100,26 @@ def get_conn_matrix( Parameters ---------- time_series : array - 2D m x n array consisting of the time-series signal for each ROI node where m = number of scans and - n = number of ROI's. + 2D m x n array consisting of the time-series signal for each ROI node + where m = number of scans and n = number of ROI's. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. dir_path : str Path to directory containing subject derivative data for given run. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's. smooth : int - Smoothing width (mm fwhm) to apply to time-series when extracting signal from ROI's. + Smoothing width (mm fwhm) to apply to time-series when extracting + signal from ROI's. dens_thresh : bool Indicates whether a target graph density is to be used as the basis for thresholding. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. ID : str A subject id or other unique identifier. roi : str @@ -149,28 +151,31 @@ def get_conn_matrix( hpass : bool High-pass filter values (Hz) to apply to node-extracted time-series. extract_strategy : str - The name of a valid function used to reduce the time-series region extraction. + The name of a valid function used to reduce the time-series region + extraction. Returns ------- conn_matrix : array Adjacency matrix stored as an m x n array of nodes and edges. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, - partcorr for partial correlation). sps type is used by default. + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial + correlation). sps type is used by default. dir_path : str Path to directory containing subject derivative data for given run. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's for tracking. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's for tracking. smooth : int - Smoothing width (mm fwhm) to apply to time-series when extracting signal from ROI's. + Smoothing width (mm fwhm) to apply to time-series when extracting + signal from ROI's. dens_thresh : bool Indicates whether a target graph density is to be used as the basis for thresholding. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. ID : str A subject id or other unique identifier. roi : str @@ -202,14 +207,17 @@ def get_conn_matrix( hpass : bool High-pass filter values (Hz) to apply to node-extracted time-series. extract_strategy : str - The name of a valid function used to reduce the time-series region extraction. + The name of a valid function used to reduce the time-series region + extraction. References ---------- - .. [1] Varoquaux, G., & Craddock, R. C. (2013). Learning and comparing functional connectomes - across subjects. NeuroImage. https://doi.org/10.1016/j.neuroimage.2013.04.007 - .. [2] Jason Laska, Manjari Narayan, 2017. skggm 0.2.7: A scikit-learn compatible package - for Gaussian and related Graphical Models. doi:10.5281/zenodo.830033 + .. [1] Varoquaux, G., & Craddock, R. C. (2013). Learning and comparing + functional connectomes across subjects. NeuroImage. + https://doi.org/10.1016/j.neuroimage.2013.04.007 + .. [2] Jason Laska, Manjari Narayan, 2017. skggm 0.2.7: + A scikit-learn compatible package for Gaussian and related Graphical + Models. doi:10.5281/zenodo.830033 """ import sys diff --git a/pynets/registration/reg_utils.py b/pynets/registration/reg_utils.py index e5148901..73df115a 100644 --- a/pynets/registration/reg_utils.py +++ b/pynets/registration/reg_utils.py @@ -294,6 +294,7 @@ def waymask2dwi_align( MNI space --> T1w --> dwi. """ import time + from nilearn.image import math_img from pynets.registration import reg_utils as regutils from nilearn.image import resample_to_img @@ -327,6 +328,10 @@ def waymask2dwi_align( time.sleep(0.5) + t_img = nib.load(waymask_in_dwi) + mask = math_img("img > 0.01", img=t_img) + mask.to_filename(waymask_in_dwi) + return waymask_in_dwi diff --git a/pynets/registration/register.py b/pynets/registration/register.py index 945fc6b9..f1ef5c76 100755 --- a/pynets/registration/register.py +++ b/pynets/registration/register.py @@ -54,8 +54,8 @@ def direct_streamline_norm( t1_aligned_mni ): """ - A Function to perform normalization of streamlines tracked in native diffusion space - to an MNI-space template. + A Function to perform normalization of streamlines tracked in native + diffusion space to an MNI-space template. Parameters ---------- @@ -66,7 +66,8 @@ def direct_streamline_norm( ap_path : str File path to the anisotropic power Nifti1Image. dir_path : str - Path to directory containing subject derivative data for a given pynets run. + Path to directory containing subject derivative data for a given + pynets run. track_type : str Tracking algorithm used (e.g. 'local' or 'particle'). target_samples : int @@ -77,8 +78,8 @@ def direct_streamline_norm( Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of brain subgraphs. node_size : int - Spherical centroid node size in the case that coordinate-based centroids - are used as ROI's for tracking. + Spherical centroid node size in the case that coordinate-based + centroids are used as ROI's for tracking. dens_thresh : bool Indicates whether a target graph density is to be used as the basis for thresholding. @@ -115,7 +116,8 @@ def direct_streamline_norm( atlas_mni : str File path to atlas parcellation Nifti1Image in T1w-warped MNI space. basedir_path : str - Path to directory to output direct-streamline normalized temp files and outputs. + Path to directory to output direct-streamline normalized temp files + and outputs. curv_thr_list : list List of integer curvature thresholds used to perform ensemble tracking. step_list : list @@ -191,8 +193,8 @@ def direct_streamline_norm( References ---------- - .. [1] Greene, C., Cieslak, M., & Grafton, S. T. (2017). Effect of different - spatial normalization approaches on tractography and structural + .. [1] Greene, C., Cieslak, M., & Grafton, S. T. (2017). Effect of + different spatial normalization approaches on tractography and structural brain networks. Network Neuroscience, 1-19. """ import gc @@ -314,21 +316,7 @@ def direct_streamline_norm( ".nii.gz", ) - # streams_warp_png = "%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s%s" % (dsn_dir, '/streamlines_mni_warp_', - # '%s' % (network + '_' if network is not - # None else ''), - # '%s' % (op.basename(roi).split('.')[0] + '_' if - # roi is not None else ''), - # conn_model, '_', target_samples, - # '%s' % ("%s%s" % - # ('_' + str(node_size), - # 'mm_') if ((node_size != 'parc') and - # (node_size is not None)) else - # '_'), - # 'curv', str(curv_thr_list).replace(', ', '_'), - # 'step', str(step_list).replace(', ', '_'), 'tracktype-', - # track_type, '_directget-', directget, '_minlength-', min_length, - # '.png') + # streams_warp_png = '/tmp/dsn.png' # SyN FA->Template [mapping, affine_map, warped_fa] = regutils.wm_syn( @@ -352,10 +340,10 @@ def direct_streamline_norm( streams_in_curr_grid = transform_streamlines( streamlines, warped_fa_affine) - # Create isocenter mapping where we anchor the origin transformation affine - # to the corner of the FOV by scaling x, y, z offsets according to a multiplicative - # van der Corput sequence with a base value equal to the voxel - # resolution + # Create isocenter mapping where we anchor the origin transformation + # affine to the corner of the FOV by scaling x, y, z offsets according + # to a multiplicative van der Corput sequence with a base value equal + # to the voxel resolution [x_mul, y_mul, z_mul] = [vdc(i, vox_size) for i in range(1, 4)] ref_grid_aff = vox_size * np.eye(4) @@ -414,8 +402,9 @@ def direct_streamline_norm( warped_fa_img.uncache() # DSN QC plotting - # plot_gen.show_template_bundles(streams_final_filt_final, atlas_mni, streams_warp_png) - # plot_gen.show_template_bundles(streamlines, fa_path, streams_warp_png) + # plot_gen.show_template_bundles(streams_final_filt_final, atlas_mni, + # streams_warp_png) plot_gen.show_template_bundles(streamlines, + # fa_path, streams_warp_png) # Create and save MNI density map nib.save( @@ -429,8 +418,8 @@ def direct_streamline_norm( density_mni, ) - # Map parcellation from native space back to MNI-space and create an 'uncertainty-union' parcellation - # with original mni-space uatlas + # Map parcellation from native space back to MNI-space and create an + # 'uncertainty-union' parcellation with original mni-space uatlas warped_uatlas = affine_map.transform_inverse( mapping.transform( np.asarray(atlas_img.dataobj).astype("int"), @@ -453,7 +442,8 @@ def direct_streamline_norm( warped_uatlas_img_res_data.astype("bool") * uatlas_mni_data.astype("bool")) os.makedirs(f"{dir_path}/parcellations", exist_ok=True) - atlas_mni = f"{dir_path}/parcellations/{op.basename(uatlas).split('.nii')[0]}_liberal.nii.gz" + atlas_mni = f"{dir_path}/parcellations/" \ + f"{op.basename(uatlas).split('.nii')[0]}_liberal.nii.gz" nib.save( nib.Nifti1Image( @@ -528,14 +518,16 @@ class DmriReg(object): References ---------- .. [1] Adluru, N., Zhang, H., Tromp, D. P. M., & Alexander, A. L. (2013). - Effects of DTI spatial normalization on white matter tract reconstructions. - Medical Imaging 2013: Image Processing. https://doi.org/10.1117/12.2007130 + Effects of DTI spatial normalization on white matter tract + reconstructions. Medical Imaging 2013: Image Processing. + https://doi.org/10.1117/12.2007130 .. [2] Greve DN, Fischl B. Accurate and robust brain image alignment using boundary-based registration. Neuroimage. 2009 Oct;48(1):63–72. doi:10.1016/j.neuroimage.2009.06.060. .. [3] Zhang Y, Brady M, Smith S. Segmentation of brain MR images through a - hidden Markov random field model and the expectation-maximization algorithm. - IEEE Trans Med Imaging. 2001 Jan;20(1):45–57. doi:10.1109/42.906424. + hidden Markov random field model and the expectation-maximization + algorithm. IEEE Trans Med Imaging. 2001 Jan;20(1):45–57. + doi:10.1109/42.906424. """ def __init__( @@ -572,32 +564,42 @@ def __init__( self.mni2t1_xfm = f"{self.reg_path_mat}{'/xfm_mni2t1.mat'}" self.mni2t1w_warp = f"{self.reg_path_warp}{'/mni2t1w_warp.nii.gz'}" self.warp_t1w2mni = f"{self.reg_path_warp}{'/t1w2mni_warp.nii.gz'}" - self.t1w2dwi = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_in_dwi.nii.gz'}" + self.t1w2dwi = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_in_dwi.nii.gz'}" self.t1_aligned_mni = ( f"{self.reg_path_img}{'/'}{self.t1w_name}{'_aligned_mni.nii.gz'}" ) - self.t1w_brain = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_brain.nii.gz'}" - self.t1w_head = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_head.nii.gz'}" + self.t1w_brain = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_brain.nii.gz'}" + self.t1w_head = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_head.nii.gz'}" self.t1w_brain_mask = ( f"{self.reg_path_img}{'/'}{self.t1w_name}{'_brain_mask.nii.gz'}" ) self.t1w_brain_mask_in_dwi = ( - f"{self.reg_path_img}{'/'}{self.t1w_name}{'_brain_mask_in_dwi.nii.gz'}" + f"{self.reg_path_img}{'/'}{self.t1w_name}" + f"{'_brain_mask_in_dwi.nii.gz'}" ) self.dwi2t1w_xfm = f"{self.reg_path_mat}{'/dwi2t1w_xfm.mat'}" self.t1w2dwi_xfm = f"{self.reg_path_mat}{'/t1w2dwi_xfm.mat'}" self.t1w2dwi_bbr_xfm = f"{self.reg_path_mat}{'/t1w2dwi_bbr_xfm.mat'}" self.dwi2t1w_bbr_xfm = f"{self.reg_path_mat}{'/dwi2t1w_bbr_xfm.mat'}" - self.t1wtissue2dwi_xfm = f"{self.reg_path_mat}{'/t1wtissue2dwi_xfm.mat'}" + self.t1wtissue2dwi_xfm = f"{self.reg_path_mat}" \ + f"{'/t1wtissue2dwi_xfm.mat'}" self.temp2dwi_xfm = ( f"{self.reg_path_mat}{'/'}{self.dwi_name}{'_xfm_temp2dwi.mat'}" ) self.map_name = f"{self.t1w_name}{'_seg'}" - self.wm_mask = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_wm.nii.gz'}" - self.wm_mask_thr = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_wm_thr.nii.gz'}" - self.wm_edge = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_wm_edge.nii.gz'}" - self.csf_mask = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_csf.nii.gz'}" - self.gm_mask = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_gm.nii.gz'}" + self.wm_mask = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_wm.nii.gz'}" + self.wm_mask_thr = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_wm_thr.nii.gz'}" + self.wm_edge = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_wm_edge.nii.gz'}" + self.csf_mask = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_csf.nii.gz'}" + self.gm_mask = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_gm.nii.gz'}" self.xfm_roi2mni_init = f"{self.reg_path_mat}{'/roi_2_mni.mat'}" self.mni_vent_loc = pkg_resources.resource_filename( "pynets", f"templates/LateralVentricles_{vox_size}.nii.gz" @@ -605,10 +607,13 @@ def __init__( self.csf_mask_dwi = ( f"{self.reg_path_img}{'/'}{self.t1w_name}{'_csf_mask_dwi.nii.gz'}" ) - self.gm_in_dwi = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_gm_in_dwi.nii.gz'}" - self.wm_in_dwi = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_wm_in_dwi.nii.gz'}" + self.gm_in_dwi = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_gm_in_dwi.nii.gz'}" + self.wm_in_dwi = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_wm_in_dwi.nii.gz'}" self.csf_mask_dwi_bin = ( - f"{self.reg_path_img}{'/'}{self.t1w_name}{'_csf_mask_dwi_bin.nii.gz'}" + f"{self.reg_path_img}{'/'}{self.t1w_name}" + f"{'_csf_mask_dwi_bin.nii.gz'}" ) self.gm_in_dwi_bin = ( f"{self.reg_path_img}{'/'}{self.t1w_name}{'_gm_in_dwi_bin.nii.gz'}" @@ -620,7 +625,8 @@ def __init__( f"{self.reg_path_img}{'/'}{self.t1w_name}{'_vent_mask_dwi.nii.gz'}" ) self.vent_csf_in_dwi = ( - f"{self.reg_path_img}{'/'}{self.t1w_name}{'_vent_csf_in_dwi.nii.gz'}" + f"{self.reg_path_img}{'/'}{self.t1w_name}" + f"{'_vent_csf_in_dwi.nii.gz'}" ) self.vent_mask_mni = f"{self.reg_path_img}{'/vent_mask_mni.nii.gz'}" self.vent_mask_t1w = f"{self.reg_path_img}{'/vent_mask_t1w.nii.gz'}" @@ -628,7 +634,7 @@ def __init__( "pynets", f"templates/{self.template_name}_{vox_size}.nii.gz" ) self.input_mni_brain = pkg_resources.resource_filename( - "pynets", f"templates/{self.template_name}_" f"brain_{vox_size}.nii.gz" + "pynets", f"templates/{self.template_name}_brain_{vox_size}.nii.gz" ) self.input_mni_mask = pkg_resources.resource_filename( "pynets", f"templates/{self.template_name}_" @@ -637,7 +643,8 @@ def __init__( "pynets", f"core/atlases/HarvardOxford-sub-prob-{vox_size}.nii.gz" ) self.wm_gm_int_in_dwi = ( - f"{self.reg_path_img}{'/'}{self.t1w_name}{'_wm_gm_int_in_dwi.nii.gz'}" + f"{self.reg_path_img}{'/'}{self.t1w_name}" + f"{'_wm_gm_int_in_dwi.nii.gz'}" ) self.wm_gm_int_in_dwi_bin = ( f"{self.reg_path_img}/{self.t1w_name}_wm_gm_int_in_dwi_bin.nii.gz" @@ -648,7 +655,8 @@ def __init__( self.corpuscallosum_mask_t1w = ( f"{self.reg_path_img}{'/CorpusCallosum_t1wmask.nii.gz'}" ) - self.corpuscallosum_dwi = f"{self.reg_path_img}{'/CorpusCallosum_dwi.nii.gz'}" + self.corpuscallosum_dwi = f"{self.reg_path_img}" \ + f"{'/CorpusCallosum_dwi.nii.gz'}" # Create empty tmp directories that do not yet exist reg_dirs = [ @@ -917,10 +925,12 @@ def t1w2dwi_align(self): def tissue2dwi_align(self): """ - A function to perform alignment of ventricle ROI's from MNI space --> dwi and CSF from T1w space --> dwi. - First generates and performs dwi space alignment of avoidance/waypoint masks for tractography. - First creates ventricle ROI. Then creates transforms from stock MNI template to dwi space. - For this to succeed, must first have called both t1w2dwi_align. + A function to perform alignment of ventricle ROI's from MNI + space --> dwi and CSF from T1w space --> dwi. First generates and + performs dwi space alignment of avoidance/waypoint masks for + tractography. First creates ventricle ROI. Then creates transforms + from stock MNI template to dwi space. For this to succeed, must first + have called both t1w2dwi_align. """ import sys import time @@ -1079,21 +1089,25 @@ def tissue2dwi_align(self): # Create ventricular CSF mask print("Creating Ventricular CSF mask...") os.system( - f"fslmaths {self.vent_mask_dwi} -kernel sphere 10 -ero -bin {self.vent_mask_dwi}" + f"fslmaths {self.vent_mask_dwi} -kernel sphere 10 -ero " + f"-bin {self.vent_mask_dwi}" ) time.sleep(1) os.system( - f"fslmaths {self.csf_mask_dwi} -add {self.vent_mask_dwi} -bin {self.vent_csf_in_dwi}" + f"fslmaths {self.csf_mask_dwi} -add {self.vent_mask_dwi} " + f"-bin {self.vent_csf_in_dwi}" ) time.sleep(1) print("Creating Corpus Callosum mask...") os.system( - f"fslmaths {self.corpuscallosum_dwi} -mas {self.wm_in_dwi_bin} -sub {self.vent_csf_in_dwi} " + f"fslmaths {self.corpuscallosum_dwi} -mas {self.wm_in_dwi_bin} " + f"-sub {self.vent_csf_in_dwi} " f"-bin {self.corpuscallosum_dwi}") time.sleep(1) # Create gm-wm interface image os.system( - f"fslmaths {self.gm_in_dwi_bin} -mul {self.wm_in_dwi_bin} -add {self.corpuscallosum_dwi} " + f"fslmaths {self.gm_in_dwi_bin} -mul {self.wm_in_dwi_bin} " + f"-add {self.corpuscallosum_dwi} " f"-mas {self.B0_mask} -bin {self.wm_gm_int_in_dwi}") time.sleep(1) return @@ -1101,15 +1115,18 @@ def tissue2dwi_align(self): class FmriReg(object): """ - A Class for Registering an atlas to a subject's MNI-aligned T1w image in native epi space. + A Class for Registering an atlas to a subject's MNI-aligned T1w image in + native epi space. + References ---------- .. [1] Brett M, Leff AP, Rorden C, Ashburner J (2001) Spatial Normalization of Brain Images with Focal Lesions Using Cost Function Masking. NeuroImage 14(2) doi:10.006/nimg.2001.0845. .. [2] Zhang Y, Brady M, Smith S. Segmentation of brain MR images through a - hidden Markov random field model and the expectation-maximization algorithm. - IEEE Trans Med Imaging. 2001 Jan;20(1):45–57. doi:10.1109/42.906424. + hidden Markov random field model and the expectation-maximization + algorithm. IEEE Trans Med Imaging. 2001 Jan;20(1):45–57. + doi:10.1109/42.906424. """ def __init__( @@ -1141,22 +1158,32 @@ def __init__( self.t1_aligned_mni = ( f"{self.reg_path_img}{'/'}{self.t1w_name}{'_aligned_mni.nii.gz'}" ) - self.t1w_brain = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_brain.nii.gz'}" - self.t1w_head = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_head.nii.gz'}" + self.t1w_brain = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_brain.nii.gz'}" + self.t1w_head = f"{self.reg_path_img}{'/'}{self.t1w_name}" \ + f"{'_head.nii.gz'}" self.t1w_brain_mask = ( - f"{self.reg_path_img}{'/'}{self.t1w_name}{'_brain_mask.nii.gz'}" + f"{self.reg_path_img}{'/'}{self.t1w_name}" + f"{'_brain_mask.nii.gz'}" ) - self.map_name = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_seg'}" - self.gm_mask = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_gm.nii.gz'}" - self.gm_mask_thr = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_gm_thr.nii.gz'}" - self.wm_mask = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_wm.nii.gz'}" - self.wm_mask_thr = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_wm_thr.nii.gz'}" - self.wm_edge = f"{self.reg_path_img}{'/'}{self.t1w_name}{'_wm_edge.nii.gz'}" + self.map_name = f"{self.reg_path_img}{'/'}" \ + f"{self.t1w_name}{'_seg'}" + self.gm_mask = f"{self.reg_path_img}{'/'}" \ + f"{self.t1w_name}{'_gm.nii.gz'}" + self.gm_mask_thr = f"{self.reg_path_img}{'/'}" \ + f"{self.t1w_name}{'_gm_thr.nii.gz'}" + self.wm_mask = f"{self.reg_path_img}{'/'}" \ + f"{self.t1w_name}{'_wm.nii.gz'}" + self.wm_mask_thr = f"{self.reg_path_img}{'/'}" \ + f"{self.t1w_name}{'_wm_thr.nii.gz'}" + self.wm_edge = f"{self.reg_path_img}{'/'}" \ + f"{self.t1w_name}{'_wm_edge.nii.gz'}" self.input_mni = pkg_resources.resource_filename( "pynets", f"templates/{self.template_name}_{vox_size}.nii.gz" ) self.input_mni_brain = pkg_resources.resource_filename( - "pynets", f"templates/{self.template_name}_" f"brain_{vox_size}.nii.gz" + "pynets", f"templates/{self.template_name}_" + f"brain_{vox_size}.nii.gz" ) self.input_mni_mask = pkg_resources.resource_filename( "pynets", f"templates/{self.template_name}_" @@ -1323,4 +1350,3 @@ def t1w2mni_align(self): self.t12mni_xfm = regutils.invert_xfm(self.mni2t1_xfm, self.t12mni_xfm) return - diff --git a/pynets/runconfig.yaml b/pynets/runconfig.yaml index 286a152d..5231c08c 100755 --- a/pynets/runconfig.yaml +++ b/pynets/runconfig.yaml @@ -30,7 +30,7 @@ tracking: sphere: - 'repulsion724' n_seeds_per_iter: # Increasing this value will decrease runtime with distributed execution, but increase runtime with serial execution. - - 1500 + - 750 max_length: - 500 tissue_classifier: # Indicates the tissue classification method to use for dMRI tractography. Options are: cmc (continuous), act (anatomically-constrained), wb (whole-brain mask), and wm (binary to white-matter only). If particle tracking is used, then this variable will be auto-set to 'cmc'. diff --git a/pynets/stats/netstats.py b/pynets/stats/netstats.py index fe01b1f7..64a049ad 100644 --- a/pynets/stats/netstats.py +++ b/pynets/stats/netstats.py @@ -113,11 +113,12 @@ def global_efficiency(G, weight="weight"): to make a completely connected graph. Since that knowlege may not exist, the scaling factor is not included. If that knowlege exists, construct the corresponding weighted graph and calculate its global_efficiency to scale - the weighted graph. Distance between nodes is calculated as the sum of weights. - If the graph is defined such that a higher weight represents a stronger connection, - distance should be represented by 1/weight. In this case, use the invert - weights function to generate a graph where the weights are set to 1/weight - and then calculate efficiency + the weighted graph. Distance between nodes is calculated as the sum of + weights. + If the graph is defined such that a higher weight represents a stronger + connection, distance should be represented by 1/weight. In this case, use + the invert weights function to generate a graph where the weights are set + to 1/weight and then calculate efficiency References ---------- @@ -244,10 +245,11 @@ def smallworldness( omega = Lr/L - C/Cl - where C and L are respectively the average clustering coefficient/transitivity - and average shortest path length of G. Lr is the average shortest path - length of an equivalent random graph and Cl is the average clustering - coefficient/transitivity of an equivalent lattice/random graph. + where C and L are respectively the average clustering + coefficient/ transitivity and average shortest path length of G. Lr is + the average shortest path length of an equivalent random graph and Cl is + the average clustering coefficient/transitivity of an equivalent + lattice/random graph. Parameters ---------- @@ -338,9 +340,11 @@ def create_communities(node_comm_aff_mat, node_num): Parameters ---------- node_comm_aff_mat : array - Community affiliation matrix produced from modularity estimation (e.g. Louvain). + Community affiliation matrix produced from modularity estimation + (e.g. Louvain). node_num : int - Number of total connected nodes in the graph used to estimate node_comm_aff_mat. + Number of total connected nodes in the graph used to estimate + node_comm_aff_mat. Returns ------- @@ -545,10 +549,10 @@ def link_communities(W, type_clustering="single"): References ---------- - .. [1] de Reus, M. A., Saenger, V. M., Kahn, R. S., & van den Heuvel, M. P. (2014). - An edge-centric perspective on the human connectome: Link communities in the - brain. Philosophical Transactions of the Royal Society B: Biological Sciences. - https://doi.org/10.1098/rstb.2013.0527 + .. [1] de Reus, M. A., Saenger, V. M., Kahn, R. S., & van den Heuvel, + M. P. (2014). An edge-centric perspective on the human connectome: + Link communities in the brain. Philosophical Transactions of the Royal + Society B: Biological Sciences. https://doi.org/10.1098/rstb.2013.0527 """ @@ -753,11 +757,11 @@ def weighted_transitivity(G): References ---------- - .. [1] Wasserman, S., and Faust, K. (1994). Social Network Analysis: Methods and - Applications. Cambridge: Cambridge University Press. - .. [2] Alain Barrat, Marc Barthelemy, Romualdo Pastor-Satorras, Alessandro Vespignani: - The architecture of complex weighted networks, Proc. Natl. Acad. Sci. - USA 101, 3747 (2004) + .. [1] Wasserman, S., and Faust, K. (1994). Social Network Analysis: + Methods and Applications. Cambridge: Cambridge University Press. + .. [2] Alain Barrat, Marc Barthelemy, Romualdo Pastor-Satorras, Alessandro + Vespignani: The architecture of complex weighted networks, Proc. Natl. + Acad. Sci. USA 101, 3747 (2004) """ @@ -828,8 +832,8 @@ def most_important(G, method="betweenness", sd=1): G : Obj NetworkX graph. method : str - Determines method for defining hubs. Valid inputs are coreness, richclub, and - eigenvector centrality. Default is coreness. + Determines method for defining hubs. Valid inputs are coreness, + richclub, and eigenvector centrality. Default is coreness. sd : int Number of standard errors as cutoff for low-importance pruning. @@ -998,13 +1002,15 @@ class CleanGraphs(object): Parameters ---------- thr : float - The value, between 0 and 1, used to threshold the graph using any variety of methods - triggered through other options. + The value, between 0 and 1, used to threshold the graph using any + variety of methods triggered through other options. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial correlation). sps type is used by default. est_path : str - File path to the thresholded graph, conn_matrix_thr, saved as a numpy array in .npy format. + File path to the thresholded graph, conn_matrix_thr, saved as a numpy + array in .npy format. prune : int Indicates whether to prune final graph of disconnected nodes/isolates. norm : int @@ -1020,9 +1026,9 @@ class CleanGraphs(object): .. [1] Qin, Tai, and Karl Rohe. "Regularized spectral clustering under the degree-corrected stochastic blockmodel." In Advances in Neural Information Processing Systems, pp. 3120-3128. 2013 - .. [2] Rohe, Karl, Tai Qin, and Bin Yu. "Co-clustering directed graphs to discover - asymmetries and directional communities." Proceedings of the National Academy - of Sciences 113.45 (2016): 12679-12684. + .. [2] Rohe, Karl, Tai Qin, and Bin Yu. "Co-clustering directed graphs to + discover asymmetries and directional communities." Proceedings of the + National Academy of Sciences 113.45 (2016): 12679-12684. """ @@ -1194,7 +1200,8 @@ def save_netmets( import os # And save results to csv out_path_neat = ( - f"{utils.create_csv_path(dir_path, est_path).split('.csv')[0]}{'_neat.csv'}" + f"{utils.create_csv_path(dir_path, est_path).split('.csv')[0]}" + f"{'_neat.csv'}" ) zipped_dict = dict(zip(metric_list_names, net_met_val_list_final)) df = pd.DataFrame.from_dict( @@ -1268,14 +1275,15 @@ def community_resolution_selection(G): resolution=resolution).values())) num_comms = len(np.unique(ci)) print( - f"{'Found '}{num_comms}{' communities at resolution: '}{resolution}{'...'}" + f"{'Found '}{num_comms}{' communities at resolution: '}" + f"{resolution}{'...'}" ) resolution = resolution + 10 tries = tries + 1 if tries > 100: print( - "\nWARNING: Louvain community detection failed. Proceeding with single community affiliation " - "vector...") + "\nWARNING: Louvain community detection failed. " + "Proceeding with single community affiliation vector...") break elif num_comms > len(G.edges()) / 10: resolution = 0.1 @@ -1288,18 +1296,20 @@ def community_resolution_selection(G): resolution=resolution).values())) num_comms = len(np.unique(ci)) print( - f"{'Found '}{num_comms}{' communities at resolution: '}{resolution}{'...'}" + f"{'Found '}{num_comms}{' communities at resolution: '}" + f"{resolution}{'...'}" ) resolution = resolution / 10 tries = tries + 1 if tries > 100: print( - "\nWARNING: Louvain community detection failed. Proceeding with single community affiliation " - "vector...") + "\nWARNING: Louvain community detection failed. " + "Proceeding with single community affiliation vector...") break else: print( - f"{'Found '}{num_comms}{' communities at resolution: '}{resolution}{'...'}" + f"{'Found '}{num_comms}{' communities at resolution: '}" + f"{resolution}{'...'}" ) return dict(zip(G.nodes(), ci)), ci, resolution, num_comms @@ -1334,14 +1344,16 @@ def get_participation(in_mat, ci, metric_list_names, net_met_val_list_final): pc_arr[j, 1] = pc_vals[j] except BaseException: print( - f"{'Participation coefficient is undefined for node '}{str(j)}{' of G'}" + f"{'Participation coefficient is undefined for node '}" + f"{str(j)}{' of G'}" ) pc_arr[j, 1] = np.nan j = j + 1 # Add mean pc_arr[num_edges, 0] = "average_participation_coefficient" nonzero_arr_partic_coef = np.delete(pc_arr[:, 1], [0]) - pc_arr[num_edges, 1] = np.nanmean(nonzero_arr_partic_coef.astype('float32'), dtype=np.float32) + pc_arr[num_edges, 1] = np.nanmean( + nonzero_arr_partic_coef.astype('float32'), dtype=np.float32) print(f"{'Mean Participation Coefficient: '}{str(pc_arr[num_edges, 1])}") for i in pc_arr[:, 0]: metric_list_names.append(i) @@ -1363,13 +1375,15 @@ def get_diversity(in_mat, ci, metric_list_names, net_met_val_list_final): dc_arr[j, 1] = dc_vals[j] except BaseException: print( - f"{'Diversity coefficient is undefined for node '}{str(j)}{' of G'}") + f"{'Diversity coefficient is undefined for node '}{str(j)}" + f"{' of G'}") dc_arr[j, 1] = np.nan j = j + 1 # Add mean dc_arr[num_edges, 0] = "average_diversity_coefficient" nonzero_arr_diversity_coef = np.delete(dc_arr[:, 1], [0]) - dc_arr[num_edges, 1] = np.nanmean(nonzero_arr_diversity_coef.astype('float32'), dtype=np.float32) + dc_arr[num_edges, 1] = np.nanmean( + nonzero_arr_diversity_coef.astype('float32'), dtype=np.float32) print(f"{'Mean Diversity Coefficient: '}{str(dc_arr[num_edges, 1])}") for i in dc_arr[:, 0]: metric_list_names.append(i) @@ -1390,12 +1404,14 @@ def get_local_efficiency(G, metric_list_names, net_met_val_list_final): try: le_arr[j, 1] = le_vals[j] except BaseException: - print(f"{'Local efficiency is undefined for node '}{str(j)}{' of G'}") + print(f"{'Local efficiency is undefined for node '}{str(j)}" + f"{' of G'}") le_arr[j, 1] = np.nan j = j + 1 le_arr[num_nodes, 0] = "average_local_efficiency_nodewise" nonzero_arr_le = np.delete(le_arr[:, 1], [0]) - le_arr[num_nodes, 1] = np.nanmean(nonzero_arr_le.astype('float32'), dtype=np.float32) + le_arr[num_nodes, 1] = np.nanmean(nonzero_arr_le.astype('float32'), + dtype=np.float32) print(f"{'Mean Local Efficiency: '}{str(le_arr[num_nodes, 1])}") for i in le_arr[:, 0]: metric_list_names.append(i) @@ -1417,12 +1433,14 @@ def get_clustering(G, metric_list_names, net_met_val_list_final): try: cl_arr[j, 1] = cl_vals[j] except BaseException: - print(f"{'Local clustering is undefined for node '}{str(j)}{' of G'}") + print(f"{'Local clustering is undefined for node '}{str(j)}" + f"{' of G'}") cl_arr[j, 1] = np.nan j = j + 1 cl_arr[num_nodes, 0] = "average_local_clustering_nodewise" nonzero_arr_cl = np.delete(cl_arr[:, 1], [0]) - cl_arr[num_nodes, 1] = np.nanmean(nonzero_arr_cl.astype('float32'), dtype=np.float32) + cl_arr[num_nodes, 1] = np.nanmean(nonzero_arr_cl.astype('float32'), + dtype=np.float32) print( f"{'Mean Local Clustering across nodes: '}{str(cl_arr[num_nodes, 1])}") for i in cl_arr[:, 0]: @@ -1446,14 +1464,17 @@ def get_degree_centrality(G, metric_list_names, net_met_val_list_final): try: dc_arr[j, 1] = dc_vals[j] except BaseException: - print(f"{'Degree centrality is undefined for node '}{str(j)}{' of G'}") + print(f"{'Degree centrality is undefined for node '}{str(j)}" + f"{' of G'}") dc_arr[j, 1] = np.nan j = j + 1 dc_arr[num_nodes, 0] = "average_degree_centrality" nonzero_arr_dc = np.delete(dc_arr[:, 1], [0]) - dc_arr[num_nodes, 1] = np.nanmean(nonzero_arr_dc.astype('float32'), dtype=np.float32) + dc_arr[num_nodes, 1] = np.nanmean(nonzero_arr_dc.astype('float32'), + dtype=np.float32) print( - f"{'Mean Degree Centrality across nodes: '}{str(dc_arr[num_nodes, 1])}") + f"{'Mean Degree Centrality across nodes: '}" + f"{str(dc_arr[num_nodes, 1])}") for i in dc_arr[:, 0]: metric_list_names.append(i) net_met_val_list_final = net_met_val_list_final + list(dc_arr[:, 1]) @@ -1479,15 +1500,18 @@ def get_betweenness_centrality( bc_arr[j, 1] = bc_vals[j] except BaseException: print( - f"{'betweennesss centrality is undefined for node '}{str(j)}{' of G'}" + f"{'betweennesss centrality is undefined for node '}" + f"{str(j)}{' of G'}" ) bc_arr[j, 1] = np.nan j = j + 1 bc_arr[num_nodes, 0] = "average_betweenness_centrality" nonzero_arr_betw_cent = np.delete(bc_arr[:, 1], [0]) - bc_arr[num_nodes, 1] = np.nanmean(nonzero_arr_betw_cent.astype('float32'), dtype=np.float32) + bc_arr[num_nodes, 1] = np.nanmean(nonzero_arr_betw_cent.astype('float32'), + dtype=np.float32) print( - f"{'Mean Betweenness Centrality across nodes: '}{str(bc_arr[num_nodes, 1])}") + f"{'Mean Betweenness Centrality across nodes: '}" + f"{str(bc_arr[num_nodes, 1])}") for i in bc_arr[:, 0]: metric_list_names.append(i) net_met_val_list_final = net_met_val_list_final + list(bc_arr[:, 1]) @@ -1510,14 +1534,17 @@ def get_eigen_centrality(G, metric_list_names, net_met_val_list_final): ec_arr[j, 1] = ec_vals[j] except BaseException: print( - f"{'Eigenvector centrality is undefined for node '}{str(j)}{' of G'}") + f"{'Eigenvector centrality is undefined for node '}" + f"{str(j)}{' of G'}") ec_arr[j, 1] = np.nan j = j + 1 ec_arr[num_nodes, 0] = "average_eigenvector_centrality" nonzero_arr_eig_cent = np.delete(ec_arr[:, 1], [0]) - ec_arr[num_nodes, 1] = np.nanmean(nonzero_arr_eig_cent.astype('float32'), dtype=np.float32) + ec_arr[num_nodes, 1] = np.nanmean(nonzero_arr_eig_cent.astype('float32'), + dtype=np.float32) print( - f"{'Mean Eigenvector Centrality across nodes: '}{str(ec_arr[num_nodes, 1])}") + f"{'Mean Eigenvector Centrality across nodes: '}" + f"{str(ec_arr[num_nodes, 1])}") for i in ec_arr[:, 0]: metric_list_names.append(i) net_met_val_list_final = net_met_val_list_final + list(ec_arr[:, 1]) @@ -1540,15 +1567,18 @@ def get_comm_centrality(G, metric_list_names, net_met_val_list_final): cc_arr[j, 1] = cc_vals[j] except BaseException: print( - f"{'Communicability centrality is undefined for node '}{str(j)}{' of G'}" + f"{'Communicability centrality is undefined for node '}" + f"{str(j)}{' of G'}" ) cc_arr[j, 1] = np.nan j = j + 1 cc_arr[num_nodes, 0] = "average_communicability_centrality" nonzero_arr_comm_cent = np.delete(cc_arr[:, 1], [0]) - cc_arr[num_nodes, 1] = np.nanmean(nonzero_arr_comm_cent.astype('float32'), dtype=np.float32) + cc_arr[num_nodes, 1] = np.nanmean(nonzero_arr_comm_cent.astype('float32'), + dtype=np.float32) print( - f"{'Mean Communicability Centrality across nodes: '}{str(cc_arr[num_nodes, 1])}" + f"{'Mean Communicability Centrality across nodes: '}" + f"{str(cc_arr[num_nodes, 1])}" ) for i in cc_arr[:, 0]: metric_list_names.append(i) @@ -1573,15 +1603,18 @@ def get_rich_club_coeff(G, metric_list_names, net_met_val_list_final): rc_arr[j, 1] = rc_vals[j] except BaseException: print( - f"{'Rich club coefficient is undefined for node '}{str(j)}{' of G'}") + f"{'Rich club coefficient is undefined for node '}" + f"{str(j)}{' of G'}") rc_arr[j, 1] = np.nan j = j + 1 # Add mean rc_arr[num_edges, 0] = "average_rich_club_coefficient" nonzero_arr_rich_club = np.delete(rc_arr[:, 1], [0]) - rc_arr[num_edges, 1] = np.nanmean(nonzero_arr_rich_club.astype('float32'), dtype=np.float32) + rc_arr[num_edges, 1] = np.nanmean(nonzero_arr_rich_club.astype('float32'), + dtype=np.float32) print( - f"{'Mean Rich Club Coefficient across edges: '}{str(rc_arr[num_edges, 1])}") + f"{'Mean Rich Club Coefficient across edges: '}" + f"{str(rc_arr[num_edges, 1])}") for i in rc_arr[:, 0]: metric_list_names.append(i) net_met_val_list_final = net_met_val_list_final + list(rc_arr[:, 1]) @@ -1606,16 +1639,18 @@ def extractnetstats( ID : str A subject id or other unique identifier. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the study of - brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. thr : float - The value, between 0 and 1, used to threshold the graph using any variety of methods - triggered through other options. + The value, between 0 and 1, used to threshold the graph using any + variety of methods triggered through other options. conn_model : str - Connectivity estimation model (e.g. corr for correlation, cov for covariance, sps for precision covariance, + Connectivity estimation model (e.g. corr for correlation, cov for + covariance, sps for precision covariance, partcorr for partial correlation). sps type is used by default. est_path : str - File path to the thresholded graph, conn_matrix_thr, saved as a numpy array in .npy format. + File path to the thresholded graph, conn_matrix_thr, saved as a numpy + array in .npy format. roi : str File path to binarized/boolean region-of-interest Nifti1Image file. prune : int @@ -1634,12 +1669,13 @@ def extractnetstats( References ---------- .. [1] Fornito, A., Zalesky, A., & Bullmore, E. T. (2016). - Fundamentals of Brain Network Analysis. In Fundamentals of Brain Network Analysis. - https://doi.org/10.1016/C2012-0-06036-X - .. [2] Aric A. Hagberg, Daniel A. Schult and Pieter J. Swart, “Exploring network structure, - dynamics, and function using NetworkX”, in Proceedings of the 7th Python in Science Conference - (SciPy2008), Gäel Varoquaux, Travis Vaught, and Jarrod Millman (Eds), (Pasadena, CA USA), - pp. 11–15, Aug 2008 + Fundamentals of Brain Network Analysis. In Fundamentals of Brain Network + Analysis. https://doi.org/10.1016/C2012-0-06036-X + .. [2] Aric A. Hagberg, Daniel A. Schult and Pieter J. Swart, + “Exploring network structure, dynamics, and function using NetworkX”, + in Proceedings of the 7th Python in Science Conference (SciPy2008), + Gäel Varoquaux, Travis Vaught, and Jarrod Millman (Eds), + (Pasadena, CA USA), pp. 11–15, Aug 2008 """ import time @@ -1766,8 +1802,8 @@ def extractnetstats( try: if ci is None: raise KeyError( - "Participation coefficient cannot be calculated for G in the absence of a " - "community affiliation vector") + "Participation coefficient cannot be calculated for G in" + " the absence of a community affiliation vector") start_time = time.time() metric_list_names, net_met_val_list_final = get_participation( in_mat, ci, metric_list_names, net_met_val_list_final @@ -1783,8 +1819,8 @@ def extractnetstats( try: if ci is None: raise KeyError( - "Diversity coefficient cannot be calculated for G in the absence of a community " - "affiliation vector") + "Diversity coefficient cannot be calculated for G in the" + " absence of a community affiliation vector") start_time = time.time() metric_list_names, net_met_val_list_final = get_diversity( in_mat, ci, metric_list_names, net_met_val_list_final @@ -1913,7 +1949,8 @@ def collect_pandas_df_make( sql_out=False, ): """ - Summarize list of pickled pandas dataframes of graph metrics unique to eacho unique combination of hyperparameters. + Summarize list of pickled pandas dataframes of graph metrics unique to + each unique combination of hyperparameters. Parameters ---------- @@ -1922,8 +1959,8 @@ def collect_pandas_df_make( ID : str A subject id or other unique identifier. network : str - Resting-state network based on Yeo-7 and Yeo-17 naming (e.g. 'Default') used to filter nodes in the - study of brain subgraphs. + Resting-state network based on Yeo-7 and Yeo-17 naming + (e.g. 'Default') used to filter nodes in the study of brain subgraphs. plot_switch : bool Activate summary plotting (histograms, central tendency, AUC, etc.) sql_out : bool @@ -1936,10 +1973,10 @@ def collect_pandas_df_make( References ---------- - .. [1] Drakesmith, M., Caeyenberghs, K., Dutt, A., Lewis, G., David, A. S., & - Jones, D. K. (2015). Overcoming the effects of false positives and threshold - bias in graph theoretical analyses of neuroimaging data. NeuroImage. - https://doi.org/10.1016/j.neuroimage.2015.05.011 + .. [1] Drakesmith, M., Caeyenberghs, K., Dutt, A., Lewis, G., David, + A. S., & Jones, D. K. (2015). Overcoming the effects of false positives + and threshold bias in graph theoretical analyses of neuroimaging data. + NeuroImage. https://doi.org/10.1016/j.neuroimage.2015.05.011 """ import sys @@ -2031,7 +2068,8 @@ def sort_thr(model_name): node_cols = [ s for s in list(df.columns) - if isinstance(s, int) or any(c.isdigit() for c in s) + if isinstance(s, int) or any(c.isdigit() for c in + s) ] if nc_collect is False: df = df.drop(node_cols, axis=1)