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/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 778817a..07decdf 100644 --- a/halfpipe2bids/main.py +++ b/halfpipe2bids/main.py @@ -1,18 +1,34 @@ from __future__ import annotations -import os +import shutil import json import pandas as pd import argparse -import logging from pathlib import Path 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 +from nilearn.connectome import ConnectivityMeasure -hp2b_log = logging.getLogger("halfpipe2bids") +hp2b_log = hp2b_logger() + + +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: @@ -41,6 +57,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 +82,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,184 +90,130 @@ 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 - ) - 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, + hp2b_log.info("Create dataset-level metadata.") + hp2b_utils.create_dataset_metadata_json( + output_dir, halfpipe_spec, path_atlas_nii + ) + 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"): + + dst = hp2b_utils.get_bids_filename(src, output_dir) + if not dst.parent.exists(): + dst.parent.mkdir(parents=True, exist_ok=True) + 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 + 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 + 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) + df_coords = pd.DataFrame( + coords, columns=["x", "y", "z"], index=atlas_label.index ) - 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 = pd.concat([atlas_label, df_coords], axis=1) + atlas_label.to_csv( + output_dir / f"{seg_meta_json.stem}.tsv", index=True, 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 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" + if args.impute_nan: + hp2b_log.info("Impute NaN with grand mean per TR.") + parcel_removal_threshold = 0.5 + 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, keep, drop = hp2b_utils.find_bad_rois( + timeseries_paths, atlas_label, parcel_removal_threshold + ) + hp2b_log.info( + "add nan imputation related information to the segmentation " + "meta data" + ) + 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" ) - 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" + 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 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.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") ) - ) - 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() - max_fd = df_conf["framewise_displacement"].max() - 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, - "MaxFramewiseDisplacement": max_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) + relmat = relmat_calculation[relmat_type].fit_transform( + [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") def main(argv: None | Sequence[str] = None) -> None: 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 49dcf8a..7cbd243 100644 --- a/halfpipe2bids/utils.py +++ b/halfpipe2bids/utils.py @@ -1,13 +1,60 @@ import os import json import pandas as pd -import logging -from nilearn.signal import clean -from nilearn import plotting import re - -hp2b_log = logging.getLogger("halfpipe2bids") -hp2b_url = "https://github.com/pbergeret12/HalfPipe2Bids/" +from halfpipe2bids import __version__ + +from halfpipe2bids.logger import hp2b_logger + +hp2b_log = hp2b_logger() +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): @@ -19,31 +66,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): @@ -62,111 +94,204 @@ 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 +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. - if df_filled.isna().any().any(): - hp2b_log.warning("Certaines valeurs n'ont pas pu être imputées.") + 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 - cleaned = clean( - df_filled.values, detrend=True, standardize="zscore_sample" + 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) + + 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 + + df_nan_prop = pd.DataFrame(per_roi_nan_counter).T / total_subjects + df_nan_prop.columns = ["proportion_missing_in_dataset"] + labels_to_drop = ( + df_nan_prop[df_nan_prop > parcel_removal_threshold] + .dropna() + .index.tolist() ) - return pd.DataFrame(cleaned, columns=df.columns, index=df.index) + 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 remove_bad_rois(dict_timeseries, label_schaefer, threshold=0.5): - # TODO: documentation - nan_counts = {label: 0 for label in label_schaefer} - total_subjects = len(dict_timeseries) +def create_dataset_metadata_json( + output_dir, halfpipe_spec, path_atlas_nii +) -> 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. + """ + # 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}") + + seg_meta = { + "File": entry + for entry in halfpipe_spec["files"] + if entry.get("suffix", False) + } - 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 + with open( + output_dir / f"seg-{seg_meta['File']['tags']['desc']}.json", "w" + ) as f: + json.dump(seg_meta, f, indent=4) - df_nan_prop = pd.DataFrame( - { - "ROI": list(nan_counts.keys()), - "proportion_nan": [ - nan_counts[label] / total_subjects for label in label_schaefer - ], - } + +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 - 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" - ) +def get_bids_filename(src, output_dir): + """ + Generates a BIDS-compliant filename based on the source file's name + and the output directory. + + 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. - return labels_to_drop + 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/. + """ + # 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] -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"] + if entities.get("sub", False): + output_dir = output_dir / f"sub-{entities['sub']}" / "func" + + if extension == ".gz": + return 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 df_coords[~df_coords.index.isin(labels_to_drop)] + return output_dir / f"{new_basename}{new_suffix_info}" -def crearte_dataset_metadata_json(output_dir) -> None: - """ - Create dataset-level metadata JSON files for BIDS. +def populate_timeseries_json(path_timeseries_json, fmriprep_dir): + """Add additional meta data for denoising metric calculation to the + existing json file. + Args: - output_dir (Path): path to the output directory where the JSON file - will be saved. + path_timeseries_json (Path): Path to the meta data file. + fmriprep_dir (Path): Associated fmriprep directory. + + Returns: + None """ - # 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, - ) + 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) - 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}", - } + sampling_freq = timeseries_meta.get("SamplingFrequency", None) - output_filename = "dataset_description.json" - output_file = output_dir / output_filename + # 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 - # Exporter le JSON - with open(output_file, "w") as f: - json.dump(json_dataset_description, f, indent=4) + # convert confound regressors + denoise_setting = timeseries_meta["Setting"] - hp2b_log.info(f"JSON exporté vers {output_dir}") + 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() + 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) diff --git a/pyproject.toml b/pyproject.toml index 119358e..b54d202 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -13,6 +13,8 @@ dependencies = [ "numpy>=2.2.6", "pandas>=2.2.3", "pip>=25.1.1", + "rich", + "tqdm>=4.67.1", ] dynamic = ["version"]