From 3881217198710f6b2b983dd065a16d33f19781ff Mon Sep 17 00:00:00 2001 From: Hao-Ting Wang Date: Wed, 2 Jul 2025 10:17:23 -0400 Subject: [PATCH 1/8] update the minimal example for set up --- CONTRIBUTING.md | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/CONTRIBUTING.md b/CONTRIBUTING.md index 68b01e4..a66fc91 100644 --- a/CONTRIBUTING.md +++ b/CONTRIBUTING.md @@ -10,8 +10,8 @@ git clone git@github.com:/HalfPipe2Bids.git Minimal example: ```bash -python3 -m venv halfpipe2bids -source halfpipe2bids/bin/activate +python3 -m venv .venv +source .venv/bin/activate ``` With `uv`: From 14a99a86742d5db93587bec138c1c6665f128b8c Mon Sep 17 00:00:00 2001 From: Hao-Ting Wang Date: Wed, 2 Jul 2025 13:21:45 -0400 Subject: [PATCH 2/8] WIP rewrite file naming --- halfpipe2bids/utils.py | 157 ++++++++++++++++++++++++++++----------- halfpipe2bids/wipmain.py | 148 ++++++++++++++++++++++++++++++++++++ 2 files changed, 260 insertions(+), 45 deletions(-) create mode 100644 halfpipe2bids/wipmain.py diff --git a/halfpipe2bids/utils.py b/halfpipe2bids/utils.py index 3d9cd79..eebb579 100644 --- a/halfpipe2bids/utils.py +++ b/halfpipe2bids/utils.py @@ -4,9 +4,58 @@ import logging from nilearn.signal import clean from nilearn import plotting +import re +from halfpipe2bids import __version__ hp2b_log = logging.getLogger("halfpipe2bids") -hp2b_url = "https://github.com/pbergeret12/HalfPipe2Bids/" +hp2b_url = "https://github.com/LAB-BRIGHT/HalfPipe2Bids" + +suffix_converter = {"matrix": "relmat", "timeseries": "timeseries"} +measure_entity_converter = { + "correlation": "PearsonCorrelation", + "covariance": "covariance", +} +regex_bids_entity = r"([a-zA-Z]*)-([^_]*)" + +dataset_description = { + "BIDSVersion": "1.9.0", + "License": None, + "Name": None, + "ReferencesAndLinks": [], + "DatasetDOI": None, + "DatasetType": "derivative", + "GeneratedBy": [ + { + "Name": "Halfpipe2Bids", + "Version": __version__, + "CodeURL": hp2b_url, + } + ], + "HowToAcknowledge": f"Please refer to our repository: {hp2b_url}", +} + +meas_meta = { + "covariates": { + "Measure": "Covariance", + "MeasureDescription": "Covariance", + "Weighted": False, + "Directed": False, + "ValidDiagonal": True, + "StorageFormat": "Full", + "NonNegative": "", + "Code": "HALFPipe", + }, + "PearsonCorrelation": { + "Measure": "Pearson correlation", + "MeasureDescription": "Pearson correlation", + "Weighted": False, + "Directed": False, + "ValidDiagonal": True, + "StorageFormat": "Full", + "NonNegative": "", + "Code": "HALFPipe", + }, +} def get_subjects(path_halfpipe_timeseries): @@ -99,57 +148,75 @@ def get_coords(volume_path, label_schaefer, labels_to_drop): return df_coords[~df_coords.index.isin(labels_to_drop)] -def crearte_dataset_metadata_json(output_dir) -> None: +def create_dataset_metadata_json(output_dir) -> None: """ Create dataset-level metadata JSON files for BIDS. Args: output_dir (Path): path to the output directory where the JSON file will be saved. """ - # export json file of common metadata for BIDS dataset - summary_path = output_dir / "meas-PearsonCorrelation_relmat.json" - with open(summary_path, "w") as f: - json.dump( - { - "Measure": "Pearson correlation", - "MeasureDescription": "Pearson correlation", - "Weighted": False, - "Directed": False, - "ValidDiagonal": True, - "StorageFormat": "Full", - "NonNegative": "", - "Code": hp2b_url, - }, - f, - indent=4, - ) + # create the dataset_description.json file + hp2b_log.info(f"Creating {output_dir / 'dataset_description.json'}") + with open(output_dir / "dataset_description.json", "w") as f: + json.dump(dataset_description, f, indent=4) + + for meas in meas_meta: + meas_path = output_dir / f"meas-{meas}_relmat.json" + with open(meas_path, "w") as f: + json.dump(meas_meta[meas], f, indent=4) + hp2b_log.info(f"Exported {meas} metadata to {meas_path}") - hp2b_log.info(f"Export terminé dans : {output_dir}") - - # Export du json de description de dataset - - json_dataset_description = { - "BIDSVersion": "1.9.0", - "License": None, - "Name": None, - "ReferencesAndLinks": [], - "DatasetDOI": None, - "DatasetType": "derivative", - "GeneratedBy": [ - { - "Name": "Halfpipe2Bids", - "Version": "0.1", - "CodeURL": hp2b_url, - } - ], - "HowToAcknowledge": f"Please refer to our repository: {hp2b_url}", - } - output_filename = "dataset_description.json" - output_file = output_dir / output_filename +def get_bids_filename(src, output_dir): + """ + Generates a BIDS-compliant filename based on the source file's name + and the output directory. - # Exporter le JSON - with open(output_file, "w") as f: - json.dump(json_dataset_description, f, indent=4) + Args: + src (Path): The source file path, which should contain BIDS + entities in its name. + output_dir (Path): The output directory where the BIDS file + will be saved. + + Returns: + Path: The BIDS-compliant file path. + + Raises: + KeyError: If required BIDS entities (e.g., 'sub', 'task', + 'atlas', 'feature') are missing from the source filename. + + Notes: + - The function extracts BIDS entities from the source filename + using a regular expression. + - It applies entity and suffix conversions according to BIDS + conventions. + - The output path is structured as: + /sub-/func/. + """ - hp2b_log.info(f"JSON exporté vers {output_dir}") + # rename files to match BIDS naming conventions + entities = re.findall(regex_bids_entity, src.stem) + entities = {entity[0]: entity[1] for entity in entities} + extension = src.suffix + suffix = src.stem.split("_")[-1] + + file_output_dir = output_dir / f"sub-{entities['sub']}" / "func" + if extension == ".gz": + return file_output_dir / src.name + + if suffix in suffix_converter: + suffix = suffix_converter[suffix] + if "desc" in entities and entities["desc"] in measure_entity_converter: + entities["desc"] = measure_entity_converter[entities["desc"]] + + # convert entities to a dictionary + new_basename = ( + f"sub-{entities['sub']}_task-{entities['task']}_" + f"seg-{entities['atlas']}_desc-{entities['feature']}_" + ) + new_suffix_info = ( + f"meas-{entities['desc']}_{suffix}{extension}" + if "desc" in entities + else f"{suffix}{extension}" + ) + return file_output_dir / f"{new_basename}{new_suffix_info}" diff --git a/halfpipe2bids/wipmain.py b/halfpipe2bids/wipmain.py new file mode 100644 index 0000000..61e0d90 --- /dev/null +++ b/halfpipe2bids/wipmain.py @@ -0,0 +1,148 @@ +from __future__ import annotations + +import shutil +import json +import pandas as pd +import argparse +import logging + +from pathlib import Path +from typing import Sequence + +from halfpipe2bids import __version__ +from halfpipe2bids import utils as hp2b_utils + +hp2b_log = logging.getLogger("halfpipe2bids") + + +timeseries_json_extra_keys = [ + "ConfoundRegressors", + "NumberOfVolumesDiscardedByMotionScrubbing", + "MeanFramewiseDisplacement", + "SamplingFrequency", + "ROICentroids", +] + + +def global_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + description=( + "Convert neuroimaging data from the HalfPipe format to the " + "standardized BIDS (Brain Imaging Data Structure) format." + ), + ) + parser.add_argument( + "halfpipe_dir", + action="store", + type=Path, + help="The directory with the HALFPipe output.", + ) + parser.add_argument( + "output_dir", + action="store", + type=Path, + help="The directory where the output files should be stored.", + ) + parser.add_argument( + "analysis_level", + help="Level of the analysis that will be performed. Only group" + " level is available.", + choices=["group"], + ) + parser.add_argument( + "--denoise-meta-data", + help="Additional metadata for denoising.", + action="store_true", + ) + parser.add_argument( + "--impute-nan", + help="Imputation and bad ROI removal.", + action="store_true", + ) + parser.add_argument( + "-v", + "--version", + action="version", + version=__version__, + ) + parser.add_argument( + "--verbosity", + help="Verbosity level.", + required=False, + choices=[0, 1, 2, 3], + default=2, + type=int, + nargs=1, + ) + return parser + + +def workflow(args: argparse.Namespace) -> None: + hp2b_log.info(vars(args)) + output_dir = args.output_dir + halfpipe_dir = args.halfpipe_dir + + # path_atlas = halfpipe_dir / "atlas" + path_derivatives = halfpipe_dir / "derivatives" + path_halfpipe_timeseries = path_derivatives / "halfpipe" + # path_fmriprep = path_derivatives / "fmriprep" + # path_label_nii = path_atlas / "atlas-Schaefer2018Combined_dseg.tsv" + path_halfpipe_spec = halfpipe_dir / "spec.json" + + if not output_dir.exists(): + output_dir.mkdir(parents=True, exist_ok=True) + + # Create dataset-level metadata + hp2b_utils.create_dataset_metadata_json(output_dir) + all_files = path_halfpipe_timeseries.glob("sub-*/**/*.*") + + # copy all files to the output directory + for src in all_files: + dst = hp2b_utils.get_bids_filename(src, output_dir) + if not dst.parent.exists(): + dst.parent.mkdir(parents=True, exist_ok=True) + hp2b_log.info(f"Renaming {src} to {dst}") + shutil.copy2(src, dst) # copy2 to preserve metadata + + # populate timeseries.json with extra information + with open(path_halfpipe_spec, "r") as f: + halfpipe_spec = json.load(f) + print(halfpipe_spec.keys()) + + +def populate_timeseries_json( + path_timeseries_json, fmriprep_dir, halfpipe_spec +): + sub = path_timeseries_json.stem.split("sub-")[-1].split("_")[0] + task = path_timeseries_json.stem.split("task-")[-1].split("_")[0] + confound_file = ( + fmriprep_dir + / f"sub-{sub}" + / "func" + / f"sub-{sub}_task-{task}_desc-confounds_timeseries.tsv" + ) + confounds = pd.read_csv(confound_file, sep="\t") + print(confounds.columns) + with open(path_timeseries_json, "r") as f: + timeseries_meta = json.load(f) + + sampling_freq = timeseries_meta.get("SamplingFrequency", None) + + # convert sampling_freq from sec to Hz + # TODO: this is an upstream issue that should be reported + if sampling_freq is not None: + sampling_freq = 1.0 / sampling_freq + + # timeseries_data.update(fmriprep_info) + # timeseries_data.update(halfpipe_info) + + with open(path_timeseries_json, "w") as f: + json.dump(timeseries_meta, f, indent=4) + + +def main(argv: None | Sequence[str] = None) -> None: + """Entry point.""" + parser = global_parser() + args = parser.parse_args(argv) + workflow(args) From 84cde27405cc7ff213db53d9600996d6c0db8352 Mon Sep 17 00:00:00 2001 From: Hao-Ting Wang Date: Wed, 2 Jul 2025 20:58:40 -0400 Subject: [PATCH 3/8] heavily refactor --- halfpipe2bids/main.py | 282 ++++++++++++++++----------------------- halfpipe2bids/utils.py | 155 +++++++++++---------- halfpipe2bids/wipmain.py | 148 -------------------- 3 files changed, 199 insertions(+), 386 deletions(-) delete mode 100644 halfpipe2bids/wipmain.py diff --git a/halfpipe2bids/main.py b/halfpipe2bids/main.py index b02cb50..e9b0fdf 100644 --- a/halfpipe2bids/main.py +++ b/halfpipe2bids/main.py @@ -1,6 +1,6 @@ from __future__ import annotations -import os +import shutil import json import pandas as pd import argparse @@ -8,6 +8,7 @@ from pathlib import Path from typing import Sequence +from nilearn.plotting import find_parcellation_cut_coords from halfpipe2bids import __version__ from halfpipe2bids import utils as hp2b_utils @@ -15,6 +16,19 @@ hp2b_log = logging.getLogger("halfpipe2bids") +def set_verbosity(verbosity: int | list[int]) -> None: + if isinstance(verbosity, list): + verbosity = verbosity[0] + if verbosity == 0: + hp2b_log.setLevel("ERROR") + elif verbosity == 1: + hp2b_log.setLevel("WARNING") + elif verbosity == 2: + hp2b_log.setLevel("INFO") + elif verbosity == 3: + hp2b_log.setLevel("DEBUG") + + def global_parser() -> argparse.ArgumentParser: parser = argparse.ArgumentParser( formatter_class=argparse.RawTextHelpFormatter, @@ -41,6 +55,16 @@ def global_parser() -> argparse.ArgumentParser: " level is available.", choices=["group"], ) + parser.add_argument( + "--denoise-metadata", + help="Add extra metadata about denoising info.", + action="store_true", + ) + parser.add_argument( + "--impute-nan", + help="Imputation and bad ROI removal.", + action="store_true", + ) parser.add_argument( "-v", "--version", @@ -56,11 +80,6 @@ def global_parser() -> argparse.ArgumentParser: type=int, nargs=1, ) - parser.add_argument( - "--NaN_Handling", - help="Enable NaN handling (imputation and bad ROI removal).", - action="store_true", - ) return parser @@ -69,182 +88,111 @@ def workflow(args: argparse.Namespace) -> None: output_dir = args.output_dir halfpipe_dir = args.halfpipe_dir - path_atlas = halfpipe_dir / "atlas" path_derivatives = halfpipe_dir / "derivatives" path_halfpipe_timeseries = path_derivatives / "halfpipe" path_fmriprep = path_derivatives / "fmriprep" - path_label_nii = path_atlas / "atlas-Schaefer2018Combined_dseg.tsv" + path_atlas_label = ( + halfpipe_dir / "atlas" / "atlas-Schaefer2018Combined_dseg.tsv" + ) + path_atlas_nii = ( + halfpipe_dir / "atlas" / "atlas-Schaefer2018Combined_dseg.nii.gz" + ) path_halfpipe_spec = halfpipe_dir / "spec.json" + set_verbosity(args.verbosity) + + with open(path_halfpipe_spec, "r") as f: + halfpipe_spec = json.load(f) if not output_dir.exists(): output_dir.mkdir(parents=True, exist_ok=True) # Create dataset-level metadata - hp2b_utils.crearte_dataset_metadata_json(output_dir) - - label_atlas = hp2b_utils.load_label_schaefer(path_label_nii) - strategy_confounds = hp2b_utils.get_strategy_confounds(path_halfpipe_spec) - subjects = hp2b_utils.get_subjects(path_halfpipe_timeseries) - - task = "task-rest" # TODO: make this dynamic - atlas_name = "schaefer400" # TODO: make this dynamic - - # --- Phase 1: Load raw time series --- - - final_timeseries = {} # final_timeseries[strategy][subject] = DataFrame - raw_data_by_subject = ( - {} - ) # raw_data_by_subject[subject][strategy] = DataFrame - - for subject in subjects: - raw_data_by_subject[subject] = {} - - for strategy in strategy_confounds: - hp_path = ( - f"{path_halfpipe_timeseries}/{subject}/func/{task}/" - f"{subject}_{task}_feature-{strategy}_atlas-{atlas_name}" - "_timeseries.tsv" - ) - if not Path(hp_path).exists(): - continue - - df = pd.read_csv(hp_path, sep="\t", header=None) - if df.shape[1] != len(label_atlas): - continue - - df.columns = label_atlas - raw_data_by_subject[subject][strategy] = df - - # --- Phase 2: Optional NaN handling (imputation and ROI filtering) --- - - if args.NaN_Handling: - for strategy in strategy_confounds: - # Build subject-wise dict for each strategy - data_for_strategy = { - subject: raw_data_by_subject[subject][strategy] - for subject in raw_data_by_subject - if strategy in raw_data_by_subject[subject] - } - - labels_to_drop = hp2b_utils.remove_bad_rois( - data_for_strategy, label_atlas + hp2b_utils.create_dataset_metadata_json( + output_dir, halfpipe_spec, path_atlas_nii + ) + all_files = path_halfpipe_timeseries.glob("sub-*/**/*.*") + + # copy all files to the output directory + for src in all_files: + dst = hp2b_utils.get_bids_filename(src, output_dir) + if not dst.parent.exists(): + dst.parent.mkdir(parents=True, exist_ok=True) + hp2b_log.info(f"Renaming {src} to {dst}") + if ".tsv" == src.suffix: + mat = pd.read_csv(src, sep="\t", header=None, na_values="nan") + mat.columns += 1 # add columns and use atlas index + mat.to_csv(dst, index=False, sep="\t", na_rep="nan") + else: + shutil.copy2(src, dst) # copy2 to preserve metadata + + if args.denoise_metadata: + # populate timeseries.json with extra information + for ts_jsons in output_dir.glob("sub-*/**/*_timeseries.json"): + hp2b_utils.populate_timeseries_json( + ts_jsons, path_fmriprep, halfpipe_spec ) - remaining_labels = [ - label for label in label_atlas if label not in labels_to_drop - ] - - for subject in data_for_strategy: - df_clean = hp2b_utils.impute_and_clean( - data_for_strategy[subject] - ) - df_clean = df_clean[remaining_labels] - - if strategy not in final_timeseries: - final_timeseries[strategy] = {} - final_timeseries[strategy][subject] = df_clean - - # Computation of remaining ROIs - coords_df = hp2b_utils.get_coords( - path_atlas / "atlas-Schaefer2018Combined_dseg.nii.gz", - label_atlas, - labels_to_drop, + seg_meta_json = list(output_dir.glob("seg-*.json"))[0] + coords = find_parcellation_cut_coords(path_atlas_nii) + atlas_label = pd.read_csv( + path_atlas_label, sep="\t", header=None, index_col=0 ) - else: - for subject in raw_data_by_subject: - for strategy in raw_data_by_subject[subject]: - if strategy not in final_timeseries: - final_timeseries[strategy] = {} - final_timeseries[strategy][subject] = raw_data_by_subject[ - subject - ][strategy] - coords_df = hp2b_utils.get_coords( - path_atlas / "atlas-Schaefer2018Combined_dseg.nii.gz", - label_atlas, - [], + atlas_label.columns = ["parcel_name"] + df_coords = pd.DataFrame( + coords, columns=["x", "y", "z"], index=atlas_label.index + ) + atlas_label = pd.concat([atlas_label, df_coords], axis=1) + atlas_label["parcel_index"] = [ + i + 1 for i in range(atlas_label.shape[0]) + ] + atlas_label.to_csv( + output_dir / f"{seg_meta_json.stem}.tsv", index=False, sep="\t" ) - # --- Phase 3: Renaming and BIDS export --- - - for subject in subjects: - for strategy in strategy_confounds: - hp2b_log.info(f"Processing {subject} | strategy: {strategy}") + if args.impute_nan: + # group file per denoising strategy + atlas_label = pd.read_csv( + path_atlas_label, sep="\t", header=None, index_col=0 + ).index.tolist() + timeseries_paths = list(output_dir.glob("sub-*/**/*_timeseries.tsv")) + dataset_nan_info = hp2b_utils.find_bad_rois( + timeseries_paths, atlas_label + ) + labels_to_drop = ( + dataset_nan_info[dataset_nan_info > 0.5].dropna().index.tolist() + ) + labels_to_keep = [ + str(label) + for label in atlas_label + if str(label) not in labels_to_drop + ] + + for p in timeseries_paths: + df = pd.read_csv(p, sep="\t", header=0, na_values="nan") + row_means = df.mean(axis=1, skipna=True) + df_imputed = df.T.fillna(row_means).T + df_imputed.loc[:, labels_to_keep].to_csv( + p, index=False, sep="\t", na_rep="nan" + ) + # TODO: recreate the functional connectivity - if subject not in final_timeseries.get(strategy, {}): - continue + seg_meta_json = list(output_dir.glob("seg-*.json"))[0] + seg_meta_tsv = list(output_dir.glob("seg-*.tsv"))[0] + seg_meta_df = pd.read_csv( + seg_meta_tsv, sep="\t", header=0, index_col="parcel_index" + ) - df_ts = final_timeseries[strategy][subject] - nroi = df_ts.shape[1] - base_name = ( - f"{subject}_{task}_seg-{atlas_name}_" - f"{nroi}_desc-denoise{strategy}" - ) - subject_output = output_dir / subject / "func" - os.makedirs(subject_output, exist_ok=True) - - # Save original ROI labels before renaming - roi_labels = df_ts.columns.tolist() - - # Save time series TSV - ts_path = subject_output / f"{base_name}_timeseries.tsv" - df_ts.columns = range(nroi) # Replace ROI names with 0...N - df_ts.to_csv(ts_path, sep="\t", index=False) - - # Save correlation matrix TSV - corr = df_ts.corr(method="pearson") - conn_path = ( - subject_output - / f"{base_name}_meas-PearsonCorrelation_relmat.tsv" - ) - corr.columns = range(nroi) - corr.to_csv(conn_path, sep="\t", index=False) - - # Load and extract metadata - json_path = ( - path_halfpipe_timeseries - / subject - / "func" - / task - / ( - f"{subject}_{task}_feature-{strategy}_atlas-" - f"{atlas_name}_timeseries.json" - ) - ) - with open(json_path) as f: - meta = json.load(f) - sampling_freq = meta.get("SamplingFrequency", None) - - # convert sampling_freq from sec to Hz - if sampling_freq is not None: - sampling_freq = 1.0 / sampling_freq - - conf_path = ( - path_fmriprep - / subject - / "func" - / f"{subject}_{task}_desc-confounds_timeseries.tsv" - ) - df_conf = pd.read_csv(conf_path, sep="\t") - mean_fd = df_conf["framewise_displacement"].mean() - scrub_vols = df_conf.filter(like="motion_outlier").shape[1] - - # ROI centroids exported as dict[label, [x, y, z]] - roi_centroids = { - label: coords_df.loc[label].tolist() - for label in roi_labels - if label in coords_df.index - } - - json_data = { - "ConfoundRegressors": strategy_confounds[strategy], - "NumberOfVolumesDiscardedByMotionScrubbing": scrub_vols, - "MeanFramewiseDisplacement": mean_fd, - "SamplingFrequency": sampling_freq, - "ROICentroids": roi_centroids, - } - - json_out_path = subject_output / f"{base_name}_timeseries.json" - with open(json_out_path, "w") as f: - json.dump(json_data, f, indent=4) + with open(seg_meta_json, "r") as f: + seg_metadata = json.load(f) + seg_metadata_exta = { + "ParcelExclusionThreashold": 0.5, + "ParcelsRemoved": labels_to_drop, + } + seg_metadata.update(seg_metadata_exta) + with open(seg_meta_json, "w") as f: + json.dump(seg_metadata, f, indent=4) + dataset_nan_info.index = seg_meta_df.index + seg_meta_df = pd.concat([seg_meta_df, dataset_nan_info], axis=1) + seg_meta_df.to_csv(seg_meta_tsv, sep="\t") def main(argv: None | Sequence[str] = None) -> None: diff --git a/halfpipe2bids/utils.py b/halfpipe2bids/utils.py index 7a87057..0a4c707 100644 --- a/halfpipe2bids/utils.py +++ b/halfpipe2bids/utils.py @@ -2,8 +2,6 @@ import json import pandas as pd import logging -from nilearn.signal import clean -from nilearn import plotting import re from halfpipe2bids import __version__ @@ -67,31 +65,16 @@ def get_subjects(path_halfpipe_timeseries): ] -def load_label_schaefer(path_label_schaefer): - # TODO: documentation and eventually remove - we what this to work with - # different atlases - return list(pd.read_csv(path_label_schaefer, sep="\t", header=None)[1]) - - -def get_strategy_confounds(spec_path): +def get_halfpipe_denoise_strategy_names(spec_path): # TODO: documentation with open(spec_path, "r") as f: data = json.load(f) - setting_to_confounds = { - s["name"]: s.get("confounds_removal", []) - for s in data.get("settings", []) - } - - strategy_confounds = {} + strategy_names = [] for feature in data.get("features", []): strategy_name = feature.get("name") - setting_name = feature.get("setting") - strategy_confounds[strategy_name] = setting_to_confounds.get( - setting_name, [] - ) - - return strategy_confounds + strategy_names.append(strategy_name) + return strategy_names def regex_to_regressor(regex_confounds, confounds_columns): @@ -110,61 +93,26 @@ def regex_to_regressor(regex_confounds, confounds_columns): return [col for col in confounds_columns if pattern.fullmatch(col)] -def impute_and_clean(df): - # TODO: documentation and what's the imputation method? - row_means = df.mean(axis=1, skipna=True) - df_filled = df.T.fillna(row_means).T - - if df_filled.isna().any().any(): - hp2b_log.warning("Certaines valeurs n'ont pas pu être imputées.") - - cleaned = clean( - df_filled.values, detrend=True, standardize="zscore_sample" - ) - return pd.DataFrame(cleaned, columns=df.columns, index=df.index) - - -def remove_bad_rois(dict_timeseries, label_schaefer, threshold=0.5): +def find_bad_rois(timeseries_paths, atlas_label): # TODO: documentation - nan_counts = {label: 0 for label in label_schaefer} - total_subjects = len(dict_timeseries) - - for df in dict_timeseries.values(): - for label in label_schaefer: - if label in df.columns and df[label].isna().all(): - nan_counts[label] += 1 - - df_nan_prop = pd.DataFrame( - { - "ROI": list(nan_counts.keys()), - "proportion_nan": [ - nan_counts[label] / total_subjects for label in label_schaefer - ], - } - ) - - labels_to_drop = df_nan_prop[df_nan_prop["proportion_nan"] > threshold][ - "ROI" - ].tolist() - - for key in dict_timeseries: - dict_timeseries[key] = dict_timeseries[key].drop( - columns=labels_to_drop, errors="ignore" - ) - - return labels_to_drop + # find out how many subject all miss the same roi + per_roi_nan_counter = {str(label): [0] for label in atlas_label} + total_subjects = len(timeseries_paths) + for p in timeseries_paths: + df = pd.read_csv(p, sep="\t", header=0, index_col=0, na_values="nan") + subject_roi_missing = (pd.isna(df).sum() / df.shape[0]) == 1 + for label in df.columns[subject_roi_missing]: + per_roi_nan_counter[label][0] += 1 -def get_coords(volume_path, label_schaefer, labels_to_drop): - # TODO: documentation - coords = plotting.find_parcellation_cut_coords(volume_path) - df_coords = pd.DataFrame( - coords, index=label_schaefer, columns=["x", "y", "z"] - ) - return df_coords[~df_coords.index.isin(labels_to_drop)] + df_nan_prop = pd.DataFrame(per_roi_nan_counter).T / total_subjects + df_nan_prop.columns = ["proportion_missing_in_dataset"] + return df_nan_prop -def create_dataset_metadata_json(output_dir) -> None: +def create_dataset_metadata_json( + output_dir, halfpipe_spec, path_atlas_nii +) -> None: """ Create dataset-level metadata JSON files for BIDS. Args: @@ -182,6 +130,17 @@ def create_dataset_metadata_json(output_dir) -> None: json.dump(meas_meta[meas], f, indent=4) hp2b_log.info(f"Exported {meas} metadata to {meas_path}") + seg_meta = { + "File": entry + for entry in halfpipe_spec["files"] + if entry.get("suffix", False) + } + + with open( + output_dir / f"seg-{seg_meta['File']['tags']['desc']}.json", "w" + ) as f: + json.dump(seg_meta, f, indent=4) + def get_bids_filename(src, output_dir): """ @@ -236,3 +195,57 @@ def get_bids_filename(src, output_dir): else f"{suffix}{extension}" ) return file_output_dir / f"{new_basename}{new_suffix_info}" + + +def populate_timeseries_json( + path_timeseries_json, fmriprep_dir, halfpipe_spec +): + """Add additional meta data for denoising metric calculation to the + existing json file. + + Args: + path_timeseries_json (Path): Path to the meta data file. + fmriprep_dir (Path): Associated fmriprep directory. + halfpipe_spec (dict): HALFPipe spec.json file. + + Returns: + None + """ + sub = path_timeseries_json.stem.split("sub-")[-1].split("_")[0] + task = path_timeseries_json.stem.split("task-")[-1].split("_")[0] + confound_file = ( + fmriprep_dir + / f"sub-{sub}" + / "func" + / f"sub-{sub}_task-{task}_desc-confounds_timeseries.tsv" + ) + confounds = pd.read_csv(confound_file, sep="\t") + extra_meta = {} + with open(path_timeseries_json, "r") as f: + timeseries_meta = json.load(f) + + sampling_freq = timeseries_meta.get("SamplingFrequency", None) + + # convert sampling_freq from sec to Hz + # TODO: this is an upstream issue that should be reported + if sampling_freq is not None: + sampling_freq = 1.0 / sampling_freq + extra_meta["SamplingFrequency"] = sampling_freq + + # convert confound regressors + denoise_setting = timeseries_meta["Setting"] + + extra_meta["ConfoundRegressors"] = regex_to_regressor( + denoise_setting["ConfoundsRemoval"], confounds.columns.tolist() + ) + extra_meta["NumberOfVolumesDiscardedByMotionScrubbing"] = len( + regex_to_regressor( + ["motion_outlier[0-9]+"], confounds.columns.tolist() + ) + ) + extra_meta["MeanFramewiseDisplacement"] = confounds[ + "framewise_displacement" + ].mean() + timeseries_meta.update(extra_meta) + with open(path_timeseries_json, "w") as f: + json.dump(timeseries_meta, f, indent=4) diff --git a/halfpipe2bids/wipmain.py b/halfpipe2bids/wipmain.py deleted file mode 100644 index 61e0d90..0000000 --- a/halfpipe2bids/wipmain.py +++ /dev/null @@ -1,148 +0,0 @@ -from __future__ import annotations - -import shutil -import json -import pandas as pd -import argparse -import logging - -from pathlib import Path -from typing import Sequence - -from halfpipe2bids import __version__ -from halfpipe2bids import utils as hp2b_utils - -hp2b_log = logging.getLogger("halfpipe2bids") - - -timeseries_json_extra_keys = [ - "ConfoundRegressors", - "NumberOfVolumesDiscardedByMotionScrubbing", - "MeanFramewiseDisplacement", - "SamplingFrequency", - "ROICentroids", -] - - -def global_parser() -> argparse.ArgumentParser: - parser = argparse.ArgumentParser( - formatter_class=argparse.RawTextHelpFormatter, - description=( - "Convert neuroimaging data from the HalfPipe format to the " - "standardized BIDS (Brain Imaging Data Structure) format." - ), - ) - parser.add_argument( - "halfpipe_dir", - action="store", - type=Path, - help="The directory with the HALFPipe output.", - ) - parser.add_argument( - "output_dir", - action="store", - type=Path, - help="The directory where the output files should be stored.", - ) - parser.add_argument( - "analysis_level", - help="Level of the analysis that will be performed. Only group" - " level is available.", - choices=["group"], - ) - parser.add_argument( - "--denoise-meta-data", - help="Additional metadata for denoising.", - action="store_true", - ) - parser.add_argument( - "--impute-nan", - help="Imputation and bad ROI removal.", - action="store_true", - ) - parser.add_argument( - "-v", - "--version", - action="version", - version=__version__, - ) - parser.add_argument( - "--verbosity", - help="Verbosity level.", - required=False, - choices=[0, 1, 2, 3], - default=2, - type=int, - nargs=1, - ) - return parser - - -def workflow(args: argparse.Namespace) -> None: - hp2b_log.info(vars(args)) - output_dir = args.output_dir - halfpipe_dir = args.halfpipe_dir - - # path_atlas = halfpipe_dir / "atlas" - path_derivatives = halfpipe_dir / "derivatives" - path_halfpipe_timeseries = path_derivatives / "halfpipe" - # path_fmriprep = path_derivatives / "fmriprep" - # path_label_nii = path_atlas / "atlas-Schaefer2018Combined_dseg.tsv" - path_halfpipe_spec = halfpipe_dir / "spec.json" - - if not output_dir.exists(): - output_dir.mkdir(parents=True, exist_ok=True) - - # Create dataset-level metadata - hp2b_utils.create_dataset_metadata_json(output_dir) - all_files = path_halfpipe_timeseries.glob("sub-*/**/*.*") - - # copy all files to the output directory - for src in all_files: - dst = hp2b_utils.get_bids_filename(src, output_dir) - if not dst.parent.exists(): - dst.parent.mkdir(parents=True, exist_ok=True) - hp2b_log.info(f"Renaming {src} to {dst}") - shutil.copy2(src, dst) # copy2 to preserve metadata - - # populate timeseries.json with extra information - with open(path_halfpipe_spec, "r") as f: - halfpipe_spec = json.load(f) - print(halfpipe_spec.keys()) - - -def populate_timeseries_json( - path_timeseries_json, fmriprep_dir, halfpipe_spec -): - sub = path_timeseries_json.stem.split("sub-")[-1].split("_")[0] - task = path_timeseries_json.stem.split("task-")[-1].split("_")[0] - confound_file = ( - fmriprep_dir - / f"sub-{sub}" - / "func" - / f"sub-{sub}_task-{task}_desc-confounds_timeseries.tsv" - ) - confounds = pd.read_csv(confound_file, sep="\t") - print(confounds.columns) - with open(path_timeseries_json, "r") as f: - timeseries_meta = json.load(f) - - sampling_freq = timeseries_meta.get("SamplingFrequency", None) - - # convert sampling_freq from sec to Hz - # TODO: this is an upstream issue that should be reported - if sampling_freq is not None: - sampling_freq = 1.0 / sampling_freq - - # timeseries_data.update(fmriprep_info) - # timeseries_data.update(halfpipe_info) - - with open(path_timeseries_json, "w") as f: - json.dump(timeseries_meta, f, indent=4) - - -def main(argv: None | Sequence[str] = None) -> None: - """Entry point.""" - parser = global_parser() - args = parser.parse_args(argv) - workflow(args) From 41814d582fd2012426d28e65431c39cb1c89b9b9 Mon Sep 17 00:00:00 2001 From: Hao-Ting Wang Date: Wed, 2 Jul 2025 21:08:27 -0400 Subject: [PATCH 4/8] fix logger --- halfpipe2bids/logger.py | 21 +++++++++++++++++++++ halfpipe2bids/main.py | 6 +++--- halfpipe2bids/utils.py | 5 +++-- pyproject.toml | 1 + 4 files changed, 28 insertions(+), 5 deletions(-) create mode 100644 halfpipe2bids/logger.py diff --git a/halfpipe2bids/logger.py b/halfpipe2bids/logger.py new file mode 100644 index 0000000..de3d78e --- /dev/null +++ b/halfpipe2bids/logger.py @@ -0,0 +1,21 @@ +"""General logger.""" + +from __future__ import annotations + +import logging + +from rich.logging import RichHandler + + +def hp2b_logger(log_level: str = "INFO") -> logging.Logger: + # FORMAT = '\n%(asctime)s - %(name)s - %(levelname)s\n\t%(message)s\n' + FORMAT = "%(message)s" + + logging.basicConfig( + level=log_level, + format=FORMAT, + datefmt="[%X]", + handlers=[RichHandler()], + ) + + return logging.getLogger("halfpipe2bids") diff --git a/halfpipe2bids/main.py b/halfpipe2bids/main.py index e9b0fdf..8fd05d8 100644 --- a/halfpipe2bids/main.py +++ b/halfpipe2bids/main.py @@ -4,7 +4,6 @@ import json import pandas as pd import argparse -import logging from pathlib import Path from typing import Sequence @@ -12,8 +11,9 @@ from halfpipe2bids import __version__ from halfpipe2bids import utils as hp2b_utils +from halfpipe2bids.logger import hp2b_logger -hp2b_log = logging.getLogger("halfpipe2bids") +hp2b_log = hp2b_logger() def set_verbosity(verbosity: int | list[int]) -> None: @@ -106,7 +106,7 @@ def workflow(args: argparse.Namespace) -> None: if not output_dir.exists(): output_dir.mkdir(parents=True, exist_ok=True) - # Create dataset-level metadata + hp2b_log.info("Create dataset-level metadata.") hp2b_utils.create_dataset_metadata_json( output_dir, halfpipe_spec, path_atlas_nii ) diff --git a/halfpipe2bids/utils.py b/halfpipe2bids/utils.py index 0a4c707..8c39897 100644 --- a/halfpipe2bids/utils.py +++ b/halfpipe2bids/utils.py @@ -1,11 +1,12 @@ import os import json import pandas as pd -import logging import re from halfpipe2bids import __version__ -hp2b_log = logging.getLogger("halfpipe2bids") +from halfpipe2bids.logger import hp2b_logger + +hp2b_log = hp2b_logger() hp2b_url = "https://github.com/LAB-BRIGHT/HalfPipe2Bids" suffix_converter = {"matrix": "relmat", "timeseries": "timeseries"} diff --git a/pyproject.toml b/pyproject.toml index 119358e..b01fa3f 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,6 +13,7 @@ dependencies = [ "numpy>=2.2.6", "pandas>=2.2.3", "pip>=25.1.1", + "rich", ] dynamic = ["version"] From d6675d066dc255e04aec2b8e837b22094aea6932 Mon Sep 17 00:00:00 2001 From: Hao-Ting Wang Date: Wed, 9 Jul 2025 11:14:53 -0400 Subject: [PATCH 5/8] draft the connectome calculation --- halfpipe2bids/main.py | 27 +++++++++++++++++++++++---- 1 file changed, 23 insertions(+), 4 deletions(-) diff --git a/halfpipe2bids/main.py b/halfpipe2bids/main.py index 8fd05d8..9359405 100644 --- a/halfpipe2bids/main.py +++ b/halfpipe2bids/main.py @@ -9,9 +9,11 @@ from typing import Sequence from nilearn.plotting import find_parcellation_cut_coords + from halfpipe2bids import __version__ from halfpipe2bids import utils as hp2b_utils from halfpipe2bids.logger import hp2b_logger +from nilearn.connectivity import ConnectivityMeasure hp2b_log = hp2b_logger() @@ -149,16 +151,19 @@ def workflow(args: argparse.Namespace) -> None: ) if args.impute_nan: - # group file per denoising strategy + parcel_removal_threshold = 0.5 atlas_label = pd.read_csv( path_atlas_label, sep="\t", header=None, index_col=0 ).index.tolist() timeseries_paths = list(output_dir.glob("sub-*/**/*_timeseries.tsv")) + # find parcels coverage stats at dataset level dataset_nan_info = hp2b_utils.find_bad_rois( timeseries_paths, atlas_label ) labels_to_drop = ( - dataset_nan_info[dataset_nan_info > 0.5].dropna().index.tolist() + dataset_nan_info[dataset_nan_info > parcel_removal_threshold] + .dropna() + .index.tolist() ) labels_to_keep = [ str(label) @@ -166,6 +171,11 @@ def workflow(args: argparse.Namespace) -> None: if str(label) not in labels_to_drop ] + # replace nan with row means (mean value of all parcels per TR) + relmat_calculation = { + "covariance": ConnectivityMeasure(kind="covariance"), + "PearsonCorrelation": ConnectivityMeasure(kind="correlation"), + } for p in timeseries_paths: df = pd.read_csv(p, sep="\t", header=0, na_values="nan") row_means = df.mean(axis=1, skipna=True) @@ -173,7 +183,16 @@ def workflow(args: argparse.Namespace) -> None: df_imputed.loc[:, labels_to_keep].to_csv( p, index=False, sep="\t", na_rep="nan" ) - # TODO: recreate the functional connectivity + # recreate the functional connectivity + for relmat_type in relmat_calculation: + dst = Path( + str(p).replace("timeseries", f"meas-{relmat_type}_relmat") + ) + relmat = relmat_calculation[relmat_type].fit_transform( + df_imputed.values + ) + df_relmat = pd.DataFrame(relmat, columns=df_imputed.columns) + df_relmat.to_csv(dst, index=False, sep="\t", na_rep="nan") seg_meta_json = list(output_dir.glob("seg-*.json"))[0] seg_meta_tsv = list(output_dir.glob("seg-*.tsv"))[0] @@ -184,7 +203,7 @@ def workflow(args: argparse.Namespace) -> None: with open(seg_meta_json, "r") as f: seg_metadata = json.load(f) seg_metadata_exta = { - "ParcelExclusionThreashold": 0.5, + "ParcelExclusionThreashold": parcel_removal_threshold, "ParcelsRemoved": labels_to_drop, } seg_metadata.update(seg_metadata_exta) From 8b9488f52c2e4fd3347e09dfaa963c70b6e0759b Mon Sep 17 00:00:00 2001 From: Hao-Ting Wang Date: Wed, 9 Jul 2025 11:17:52 -0400 Subject: [PATCH 6/8] apply changes from pr#32 --- halfpipe2bids/_oldmain.py | 254 ++++++++++++++++++++++++++++++++++++++ halfpipe2bids/utils.py | 3 + 2 files changed, 257 insertions(+) create mode 100644 halfpipe2bids/_oldmain.py diff --git a/halfpipe2bids/_oldmain.py b/halfpipe2bids/_oldmain.py new file mode 100644 index 0000000..b02cb50 --- /dev/null +++ b/halfpipe2bids/_oldmain.py @@ -0,0 +1,254 @@ +from __future__ import annotations + +import os +import json +import pandas as pd +import argparse +import logging + +from pathlib import Path +from typing import Sequence + +from halfpipe2bids import __version__ +from halfpipe2bids import utils as hp2b_utils + +hp2b_log = logging.getLogger("halfpipe2bids") + + +def global_parser() -> argparse.ArgumentParser: + parser = argparse.ArgumentParser( + formatter_class=argparse.RawTextHelpFormatter, + description=( + "Convert neuroimaging data from the HalfPipe format to the " + "standardized BIDS (Brain Imaging Data Structure) format." + ), + ) + parser.add_argument( + "halfpipe_dir", + action="store", + type=Path, + help="The directory with the HALFPipe output.", + ) + parser.add_argument( + "output_dir", + action="store", + type=Path, + help="The directory where the output files should be stored.", + ) + parser.add_argument( + "analysis_level", + help="Level of the analysis that will be performed. Only group" + " level is available.", + choices=["group"], + ) + parser.add_argument( + "-v", + "--version", + action="version", + version=__version__, + ) + parser.add_argument( + "--verbosity", + help="Verbosity level.", + required=False, + choices=[0, 1, 2, 3], + default=2, + type=int, + nargs=1, + ) + parser.add_argument( + "--NaN_Handling", + help="Enable NaN handling (imputation and bad ROI removal).", + action="store_true", + ) + return parser + + +def workflow(args: argparse.Namespace) -> None: + hp2b_log.info(vars(args)) + output_dir = args.output_dir + halfpipe_dir = args.halfpipe_dir + + path_atlas = halfpipe_dir / "atlas" + path_derivatives = halfpipe_dir / "derivatives" + path_halfpipe_timeseries = path_derivatives / "halfpipe" + path_fmriprep = path_derivatives / "fmriprep" + path_label_nii = path_atlas / "atlas-Schaefer2018Combined_dseg.tsv" + path_halfpipe_spec = halfpipe_dir / "spec.json" + + if not output_dir.exists(): + output_dir.mkdir(parents=True, exist_ok=True) + + # Create dataset-level metadata + hp2b_utils.crearte_dataset_metadata_json(output_dir) + + label_atlas = hp2b_utils.load_label_schaefer(path_label_nii) + strategy_confounds = hp2b_utils.get_strategy_confounds(path_halfpipe_spec) + subjects = hp2b_utils.get_subjects(path_halfpipe_timeseries) + + task = "task-rest" # TODO: make this dynamic + atlas_name = "schaefer400" # TODO: make this dynamic + + # --- Phase 1: Load raw time series --- + + final_timeseries = {} # final_timeseries[strategy][subject] = DataFrame + raw_data_by_subject = ( + {} + ) # raw_data_by_subject[subject][strategy] = DataFrame + + for subject in subjects: + raw_data_by_subject[subject] = {} + + for strategy in strategy_confounds: + hp_path = ( + f"{path_halfpipe_timeseries}/{subject}/func/{task}/" + f"{subject}_{task}_feature-{strategy}_atlas-{atlas_name}" + "_timeseries.tsv" + ) + if not Path(hp_path).exists(): + continue + + df = pd.read_csv(hp_path, sep="\t", header=None) + if df.shape[1] != len(label_atlas): + continue + + df.columns = label_atlas + raw_data_by_subject[subject][strategy] = df + + # --- Phase 2: Optional NaN handling (imputation and ROI filtering) --- + + if args.NaN_Handling: + for strategy in strategy_confounds: + # Build subject-wise dict for each strategy + data_for_strategy = { + subject: raw_data_by_subject[subject][strategy] + for subject in raw_data_by_subject + if strategy in raw_data_by_subject[subject] + } + + labels_to_drop = hp2b_utils.remove_bad_rois( + data_for_strategy, label_atlas + ) + remaining_labels = [ + label for label in label_atlas if label not in labels_to_drop + ] + + for subject in data_for_strategy: + df_clean = hp2b_utils.impute_and_clean( + data_for_strategy[subject] + ) + df_clean = df_clean[remaining_labels] + + if strategy not in final_timeseries: + final_timeseries[strategy] = {} + final_timeseries[strategy][subject] = df_clean + + # Computation of remaining ROIs + coords_df = hp2b_utils.get_coords( + path_atlas / "atlas-Schaefer2018Combined_dseg.nii.gz", + label_atlas, + labels_to_drop, + ) + else: + for subject in raw_data_by_subject: + for strategy in raw_data_by_subject[subject]: + if strategy not in final_timeseries: + final_timeseries[strategy] = {} + final_timeseries[strategy][subject] = raw_data_by_subject[ + subject + ][strategy] + coords_df = hp2b_utils.get_coords( + path_atlas / "atlas-Schaefer2018Combined_dseg.nii.gz", + label_atlas, + [], + ) + + # --- Phase 3: Renaming and BIDS export --- + + for subject in subjects: + for strategy in strategy_confounds: + hp2b_log.info(f"Processing {subject} | strategy: {strategy}") + + if subject not in final_timeseries.get(strategy, {}): + continue + + df_ts = final_timeseries[strategy][subject] + nroi = df_ts.shape[1] + base_name = ( + f"{subject}_{task}_seg-{atlas_name}_" + f"{nroi}_desc-denoise{strategy}" + ) + subject_output = output_dir / subject / "func" + os.makedirs(subject_output, exist_ok=True) + + # Save original ROI labels before renaming + roi_labels = df_ts.columns.tolist() + + # Save time series TSV + ts_path = subject_output / f"{base_name}_timeseries.tsv" + df_ts.columns = range(nroi) # Replace ROI names with 0...N + df_ts.to_csv(ts_path, sep="\t", index=False) + + # Save correlation matrix TSV + corr = df_ts.corr(method="pearson") + conn_path = ( + subject_output + / f"{base_name}_meas-PearsonCorrelation_relmat.tsv" + ) + corr.columns = range(nroi) + corr.to_csv(conn_path, sep="\t", index=False) + + # Load and extract metadata + json_path = ( + path_halfpipe_timeseries + / subject + / "func" + / task + / ( + f"{subject}_{task}_feature-{strategy}_atlas-" + f"{atlas_name}_timeseries.json" + ) + ) + with open(json_path) as f: + meta = json.load(f) + sampling_freq = meta.get("SamplingFrequency", None) + + # convert sampling_freq from sec to Hz + if sampling_freq is not None: + sampling_freq = 1.0 / sampling_freq + + conf_path = ( + path_fmriprep + / subject + / "func" + / f"{subject}_{task}_desc-confounds_timeseries.tsv" + ) + df_conf = pd.read_csv(conf_path, sep="\t") + mean_fd = df_conf["framewise_displacement"].mean() + scrub_vols = df_conf.filter(like="motion_outlier").shape[1] + + # ROI centroids exported as dict[label, [x, y, z]] + roi_centroids = { + label: coords_df.loc[label].tolist() + for label in roi_labels + if label in coords_df.index + } + + json_data = { + "ConfoundRegressors": strategy_confounds[strategy], + "NumberOfVolumesDiscardedByMotionScrubbing": scrub_vols, + "MeanFramewiseDisplacement": mean_fd, + "SamplingFrequency": sampling_freq, + "ROICentroids": roi_centroids, + } + + json_out_path = subject_output / f"{base_name}_timeseries.json" + with open(json_out_path, "w") as f: + json.dump(json_data, f, indent=4) + + +def main(argv: None | Sequence[str] = None) -> None: + """Entry point.""" + parser = global_parser() + args = parser.parse_args(argv) + workflow(args) diff --git a/halfpipe2bids/utils.py b/halfpipe2bids/utils.py index 8c39897..740ae51 100644 --- a/halfpipe2bids/utils.py +++ b/halfpipe2bids/utils.py @@ -247,6 +247,9 @@ def populate_timeseries_json( extra_meta["MeanFramewiseDisplacement"] = confounds[ "framewise_displacement" ].mean() + extra_meta["MaxFramewiseDisplacement"] = ( + confounds["framewise_displacement"].max(), + ) timeseries_meta.update(extra_meta) with open(path_timeseries_json, "w") as f: json.dump(timeseries_meta, f, indent=4) From c0dc3b1eb5c297f3e6d732293f74182b1f2ea57b Mon Sep 17 00:00:00 2001 From: Hao-Ting Wang Date: Wed, 9 Jul 2025 11:22:31 -0400 Subject: [PATCH 7/8] wrong module name --- halfpipe2bids/main.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/halfpipe2bids/main.py b/halfpipe2bids/main.py index 9359405..d137068 100644 --- a/halfpipe2bids/main.py +++ b/halfpipe2bids/main.py @@ -13,7 +13,7 @@ from halfpipe2bids import __version__ from halfpipe2bids import utils as hp2b_utils from halfpipe2bids.logger import hp2b_logger -from nilearn.connectivity import ConnectivityMeasure +from nilearn.connectome import ConnectivityMeasure hp2b_log = hp2b_logger() @@ -200,6 +200,7 @@ def workflow(args: argparse.Namespace) -> None: seg_meta_tsv, sep="\t", header=0, index_col="parcel_index" ) + # add nan imputation related information to the segmentation meta data with open(seg_meta_json, "r") as f: seg_metadata = json.load(f) seg_metadata_exta = { From cbabba166dee18e5fe737ee41fbcd75636878c5c Mon Sep 17 00:00:00 2001 From: Hao-Ting Wang Date: Fri, 11 Jul 2025 17:52:38 -0400 Subject: [PATCH 8/8] basic tests --- halfpipe2bids/main.py | 123 ++++++++++++++++---------------- halfpipe2bids/tests/test_cli.py | 52 +++++++++----- halfpipe2bids/utils.py | 64 ++++++++++++++--- pyproject.toml | 1 + 4 files changed, 150 insertions(+), 90 deletions(-) diff --git a/halfpipe2bids/main.py b/halfpipe2bids/main.py index d137068..07decdf 100644 --- a/halfpipe2bids/main.py +++ b/halfpipe2bids/main.py @@ -9,7 +9,7 @@ from typing import Sequence from nilearn.plotting import find_parcellation_cut_coords - +from tqdm import tqdm from halfpipe2bids import __version__ from halfpipe2bids import utils as hp2b_utils from halfpipe2bids.logger import hp2b_logger @@ -112,108 +112,109 @@ def workflow(args: argparse.Namespace) -> None: hp2b_utils.create_dataset_metadata_json( output_dir, halfpipe_spec, path_atlas_nii ) - all_files = path_halfpipe_timeseries.glob("sub-*/**/*.*") + all_files = path_halfpipe_timeseries.glob("sub-*/**/sub-*.*") + + hp2b_log.info(f"Copy all files to the output directory: {output_dir}") + for src in tqdm(all_files, desc="Renaming files"): - # copy all files to the output directory - for src in all_files: dst = hp2b_utils.get_bids_filename(src, output_dir) if not dst.parent.exists(): dst.parent.mkdir(parents=True, exist_ok=True) - hp2b_log.info(f"Renaming {src} to {dst}") - if ".tsv" == src.suffix: + hp2b_log.debug(f"Renaming {src} to {dst}") + if ".tsv" == src.suffix: # add columns and use atlas index mat = pd.read_csv(src, sep="\t", header=None, na_values="nan") - mat.columns += 1 # add columns and use atlas index + mat.columns += 1 mat.to_csv(dst, index=False, sep="\t", na_rep="nan") else: shutil.copy2(src, dst) # copy2 to preserve metadata if args.denoise_metadata: # populate timeseries.json with extra information - for ts_jsons in output_dir.glob("sub-*/**/*_timeseries.json"): - hp2b_utils.populate_timeseries_json( - ts_jsons, path_fmriprep, halfpipe_spec - ) seg_meta_json = list(output_dir.glob("seg-*.json"))[0] + all_meta_json = output_dir.glob("sub-*/**/*_timeseries.json") + for ts_jsons in all_meta_json: + hp2b_utils.populate_timeseries_json(ts_jsons, path_fmriprep) + + atlas_label = hp2b_utils.load_atlas_info_tsv(path_atlas_label) coords = find_parcellation_cut_coords(path_atlas_nii) - atlas_label = pd.read_csv( - path_atlas_label, sep="\t", header=None, index_col=0 - ) - atlas_label.columns = ["parcel_name"] df_coords = pd.DataFrame( coords, columns=["x", "y", "z"], index=atlas_label.index ) atlas_label = pd.concat([atlas_label, df_coords], axis=1) - atlas_label["parcel_index"] = [ - i + 1 for i in range(atlas_label.shape[0]) - ] atlas_label.to_csv( - output_dir / f"{seg_meta_json.stem}.tsv", index=False, sep="\t" + output_dir / f"{seg_meta_json.stem}.tsv", index=True, sep="\t" ) if args.impute_nan: + hp2b_log.info("Impute NaN with grand mean per TR.") parcel_removal_threshold = 0.5 - atlas_label = pd.read_csv( - path_atlas_label, sep="\t", header=None, index_col=0 - ).index.tolist() + seg_meta_df = hp2b_utils.load_atlas_info_tsv(path_atlas_label) + atlas_label = seg_meta_df.index.tolist() timeseries_paths = list(output_dir.glob("sub-*/**/*_timeseries.tsv")) # find parcels coverage stats at dataset level - dataset_nan_info = hp2b_utils.find_bad_rois( - timeseries_paths, atlas_label + dataset_nan_info, keep, drop = hp2b_utils.find_bad_rois( + timeseries_paths, atlas_label, parcel_removal_threshold ) - labels_to_drop = ( - dataset_nan_info[dataset_nan_info > parcel_removal_threshold] - .dropna() - .index.tolist() + hp2b_log.info( + "add nan imputation related information to the segmentation " + "meta data" ) - labels_to_keep = [ - str(label) - for label in atlas_label - if str(label) not in labels_to_drop - ] + seg_meta_json = list(output_dir.glob("seg-*.json"))[0] + seg_meta_tsv = output_dir / f"{seg_meta_json.stem}.tsv" + if seg_meta_tsv.exists(): + seg_meta_df = pd.read_csv( + seg_meta_tsv, sep="\t", header=0, index_col="parcel_index" + ) + else: + seg_meta_df = hp2b_utils.load_atlas_info_tsv(path_atlas_label) + + with open(seg_meta_json, "r") as f: + seg_metadata = json.load(f) + seg_metadata_exta = { + "ParcelExclusionThreashold": parcel_removal_threshold, + "ParcelsRemoved": drop, + } + seg_metadata.update(seg_metadata_exta) + with open(seg_meta_json, "w") as f: + json.dump(seg_metadata, f, indent=4) + dataset_nan_info.index = seg_meta_df.index + seg_meta_df = pd.concat([seg_meta_df, dataset_nan_info], axis=1) + seg_meta_df.to_csv(seg_meta_tsv, sep="\t") + + hp2b_log.info( + f"Dropping {len(seg_metadata_exta['ParcelsRemoved'])} " + f"ROIs due to {parcel_removal_threshold*100}% of the " + "subject have no signal these regions." + ) # replace nan with row means (mean value of all parcels per TR) relmat_calculation = { "covariance": ConnectivityMeasure(kind="covariance"), "PearsonCorrelation": ConnectivityMeasure(kind="correlation"), } - for p in timeseries_paths: - df = pd.read_csv(p, sep="\t", header=0, na_values="nan") - row_means = df.mean(axis=1, skipna=True) + for p in tqdm( + timeseries_paths, + desc="Imputing NaN and recalculate functional connectomes", + ): + df = pd.read_csv(p, sep="\t", header=0, na_values="nan").loc[ + :, keep + ] + row_means = df.mean(axis=1, skipna=True) # global mean per TR df_imputed = df.T.fillna(row_means).T - df_imputed.loc[:, labels_to_keep].to_csv( - p, index=False, sep="\t", na_rep="nan" - ) + df_imputed.to_csv(p, index=False, sep="\t", na_rep="nan") + hp2b_log.debug(df_imputed.shape) + hp2b_log.debug(p) # recreate the functional connectivity for relmat_type in relmat_calculation: dst = Path( str(p).replace("timeseries", f"meas-{relmat_type}_relmat") ) relmat = relmat_calculation[relmat_type].fit_transform( - df_imputed.values - ) + [df_imputed.values] + )[0] df_relmat = pd.DataFrame(relmat, columns=df_imputed.columns) df_relmat.to_csv(dst, index=False, sep="\t", na_rep="nan") - seg_meta_json = list(output_dir.glob("seg-*.json"))[0] - seg_meta_tsv = list(output_dir.glob("seg-*.tsv"))[0] - seg_meta_df = pd.read_csv( - seg_meta_tsv, sep="\t", header=0, index_col="parcel_index" - ) - - # add nan imputation related information to the segmentation meta data - with open(seg_meta_json, "r") as f: - seg_metadata = json.load(f) - seg_metadata_exta = { - "ParcelExclusionThreashold": parcel_removal_threshold, - "ParcelsRemoved": labels_to_drop, - } - seg_metadata.update(seg_metadata_exta) - with open(seg_meta_json, "w") as f: - json.dump(seg_metadata, f, indent=4) - dataset_nan_info.index = seg_meta_df.index - seg_meta_df = pd.concat([seg_meta_df, dataset_nan_info], axis=1) - seg_meta_df.to_csv(seg_meta_tsv, sep="\t") - def main(argv: None | Sequence[str] = None) -> None: """Entry point.""" diff --git a/halfpipe2bids/tests/test_cli.py b/halfpipe2bids/tests/test_cli.py index e26cd94..bfa53b0 100644 --- a/halfpipe2bids/tests/test_cli.py +++ b/halfpipe2bids/tests/test_cli.py @@ -39,37 +39,53 @@ def test_smoke(tmp_path, caplog): / "tests/data/dataset-ds000030_halfpipe1.2.3dev" ) output_dir = tmp_path / "output" - - main( - [ - str(halfpipe_dir), - str(output_dir), - "group", - ] - ) + cmd = [ + str(halfpipe_dir), + str(output_dir), + "group", + ] + main(cmd) output_folder = output_dir / "sub-10159/func" - base = "sub-10159_task-rest_seg-schaefer400_434" - ts_base = base + "_desc-denoisecorrMatrix1" + base = "sub-10159_task-rest_seg-schaefer400" + ts_base = base + "_desc-corrMatrix1" relmat_file = output_folder / ( ts_base + "_meas-PearsonCorrelation_relmat.tsv" ) # checking if relmat file exists - if not relmat_file.exists(): - raise FileNotFoundError( - f"Expected file not found: {relmat_file}\nAvailable files:\n" - + "\n".join(str(p) for p in output_folder.glob("*")) - ) - + assert relmat_file.exists() relmat = pd.read_csv(relmat_file, sep="\t") # This is the number of ROI (columns) I got from the supposedly original file - # TODO: when the --impute-nans option is added, this test should pass assert relmat.shape[1] == 434 json_file = output_folder / (ts_base + "_timeseries.json") assert json_file.exists() with open(json_file, "r") as f: content = json.load(f) # the unit is Hz, for TR= 2s, the sampling frequency is 0.5 Hz - assert content.get("SamplingFrequency") == 0.5 + # however, when no flags are passed, this is just copying the original + # file, hence ,mistake remains. + assert content.get("SamplingFrequency") == 2 # TODO: when the --impute-nans option is added, create a test for # the relmat with NaNs replaced by grand mean + + main(cmd + ["--denoise-meta"]) + assert json_file.exists() + with open(json_file, "r") as f: + content = json.load(f) + # the unit is Hz, for TR= 2s, the sampling frequency is 0.5 Hz + assert content.get("SamplingFrequency") == 0.5 + relmat = pd.read_csv(relmat_file, sep="\t") + # This is the number of ROI (columns) I got from the supposedly original file + assert relmat.shape[1] == 434 # the content of the file untouched + + main(cmd + ["--impute-nan"]) + assert json_file.exists() + with open(json_file, "r") as f: + content = json.load(f) + # the unit is Hz, for TR= 2s, the sampling frequency is 0.5 Hz + # however, when no flags are passed, this is just copying the original + # file, hence ,mistake remains. + assert content.get("SamplingFrequency") == 2 + relmat = pd.read_csv(relmat_file, sep="\t") + # This is the number of ROI (columns) I got from the supposedly original file + assert relmat.shape[1] == 417 # ROI with too many subjects missing removed diff --git a/halfpipe2bids/utils.py b/halfpipe2bids/utils.py index 740ae51..7cbd243 100644 --- a/halfpipe2bids/utils.py +++ b/halfpipe2bids/utils.py @@ -94,9 +94,24 @@ def regex_to_regressor(regex_confounds, confounds_columns): return [col for col in confounds_columns if pattern.fullmatch(col)] -def find_bad_rois(timeseries_paths, atlas_label): - # TODO: documentation - # find out how many subject all miss the same roi +def find_bad_rois(timeseries_paths, atlas_label, parcel_removal_threshold=0.5): + """ + Find out how many subject miss the same roi report in proportion of the + dataset. + + Args: + timeseries_paths (list[Path]): Path to all time series data. + atlas_label (List): Parcel index (starting from 1). + parcel_removal_threshold (float): proportion of the dataset. + 1.0 = all subjects in the dataset miss a given parcel. + 0.0 = all subjects in the dataset has a given parcel. + Default: 0.5 + + Returns: + pandas.DataFrame: proportion of the dataset with nan per parcel. + List: labels to keep. + List: labels to drop. + """ per_roi_nan_counter = {str(label): [0] for label in atlas_label} total_subjects = len(timeseries_paths) @@ -108,7 +123,15 @@ def find_bad_rois(timeseries_paths, atlas_label): df_nan_prop = pd.DataFrame(per_roi_nan_counter).T / total_subjects df_nan_prop.columns = ["proportion_missing_in_dataset"] - return df_nan_prop + labels_to_drop = ( + df_nan_prop[df_nan_prop > parcel_removal_threshold] + .dropna() + .index.tolist() + ) + labels_to_keep = [ + str(label) for label in atlas_label if str(label) not in labels_to_drop + ] + return df_nan_prop, labels_to_keep, labels_to_drop def create_dataset_metadata_json( @@ -143,6 +166,26 @@ def create_dataset_metadata_json( json.dump(seg_meta, f, indent=4) +def load_atlas_info_tsv(path_atlas_label): + """ + Load original atlas parcel label and index tsv from halfpipe. + The first column is the index, second the parcel label. + There's no header in the file. + + Args: + path_atlas_label (Path): Path to the file. + + Returns: + pandas.DataFrame: + """ + atlas_label = pd.read_csv( + path_atlas_label, sep="\t", header=None, index_col=0 + ) + atlas_label.columns = ["parcel_name"] + atlas_label.index.name = "parcel_index" + return atlas_label + + def get_bids_filename(src, output_dir): """ Generates a BIDS-compliant filename based on the source file's name @@ -176,9 +219,11 @@ def get_bids_filename(src, output_dir): extension = src.suffix suffix = src.stem.split("_")[-1] - file_output_dir = output_dir / f"sub-{entities['sub']}" / "func" + if entities.get("sub", False): + output_dir = output_dir / f"sub-{entities['sub']}" / "func" + if extension == ".gz": - return file_output_dir / src.name + return output_dir / src.name if suffix in suffix_converter: suffix = suffix_converter[suffix] @@ -195,19 +240,16 @@ def get_bids_filename(src, output_dir): if "desc" in entities else f"{suffix}{extension}" ) - return file_output_dir / f"{new_basename}{new_suffix_info}" + return output_dir / f"{new_basename}{new_suffix_info}" -def populate_timeseries_json( - path_timeseries_json, fmriprep_dir, halfpipe_spec -): +def populate_timeseries_json(path_timeseries_json, fmriprep_dir): """Add additional meta data for denoising metric calculation to the existing json file. Args: path_timeseries_json (Path): Path to the meta data file. fmriprep_dir (Path): Associated fmriprep directory. - halfpipe_spec (dict): HALFPipe spec.json file. Returns: None diff --git a/pyproject.toml b/pyproject.toml index b01fa3f..b54d202 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -14,6 +14,7 @@ dependencies = [ "pandas>=2.2.3", "pip>=25.1.1", "rich", + "tqdm>=4.67.1", ] dynamic = ["version"]