diff --git a/code/count_frames.py b/code/count_frames.py
new file mode 100644
index 0000000..6ecca8a
--- /dev/null
+++ b/code/count_frames.py
@@ -0,0 +1,80 @@
+import pandas as pd
+import argparse
+import h5py
+
+from pathlib import Path
+
+
+def count_frames(hdf5_object, target_suffix="atlas-MIST_desc-64_timeseries", path=""):
+ data = []
+ for name in hdf5_object:
+ item_path = f"{path}/{name}" if path else name
+ if isinstance(hdf5_object[name], h5py.Dataset) and name.endswith(target_suffix):
+ # Collect dataset name and the first dimension of its shape
+ data.append(
+ {"file": name, "total_frames": hdf5_object[name].shape[0]}
+ ) # n_frames
+ elif isinstance(hdf5_object[name], h5py.Group):
+ # Recursive call to navigate through groups and collect data
+ data.extend(count_frames(hdf5_object[name], target_suffix, item_path))
+ return data
+
+
+if __name__ == "__main__":
+ parser = argparse.ArgumentParser(
+ description="Create a .tsv listing number of frames remaining for speicific atlas/strategy"
+ )
+ parser.add_argument(
+ "input_p",
+ type=Path,
+ help="Path to directory containing subject directories [participant] or .h5 file of connectomes [group]",
+ )
+ parser.add_argument("output_p", type=Path, help="Output path")
+
+ parser.add_argument(
+ "level",
+ choices=["group", "participant"],
+ help="Analysis level: 'group' or 'participant'",
+ )
+
+ parser.add_argument("dataset", type=str, help="Name of the dataset")
+
+ args = parser.parse_args()
+
+ data = []
+ if args.level == "participant":
+ for sub_directory in args.input_p.iterdir():
+ if sub_directory.is_dir() and sub_directory.name != "working_directory":
+ hdf5_file_p = (
+ sub_directory
+ / f"sub-{sub_directory.name}_atlas-MIST_desc-simple.h5" # atlas-MIST_desc-scrubbing.5+gsr.h5"
+ )
+ if hdf5_file_p.exists():
+ with h5py.File(hdf5_file_p, "r") as file:
+ data.extend(count_frames(file))
+ else:
+ print(f"File {hdf5_file_p} does not exist.")
+ elif args.level == "group":
+ hdf5_file_p = (
+ args.input_p / "atlas-MIST_desc-simple.h5"
+ ) # "atlas-MIST_desc-scrubbing.5+gsr.h5"
+ with h5py.File(hdf5_file_p, "r") as file:
+ data = count_frames(file)
+
+ # Convert the list of dictionaries into a DataFrame
+ df = pd.DataFrame(data)
+
+ # Extract file details. Session and run are optional
+ pattern = (
+ r"sub-([^\s/_]+)(?:_ses-([a-zA-Z0-9]+))?_task-([a-zA-Z0-9]+)(?:_run-(\d+))?"
+ )
+ df[["participant_id", "ses", "task", "run"]] = df["file"].str.extract(pattern)
+ df["participant_id"] = df["participant_id"].astype(str)
+
+ # Drop the file column and add dataset as a column
+ df.drop(columns=["file"], inplace=True)
+ df["dataset"] = args.dataset
+
+ # Save output
+ # df.to_csv(args.output_p / f"{args.dataset}_frames.tsv", sep="\t", index=False)
+ df.to_csv(args.output_p / f"{args.dataset}_frames_total.tsv", sep="\t", index=False)
diff --git a/code/filter_qc.py b/code/filter_qc.py
new file mode 100644
index 0000000..b82cf1a
--- /dev/null
+++ b/code/filter_qc.py
@@ -0,0 +1,161 @@
+"""Load phenotypic and QC data for datasets, summarise QC pass rate, and output a .tsv with passed QC scans in order to obtain connectomes.
+Optionally, pass specific diagnoses.
+
+To run, needs:
+- TSV of pheno and qc data for each specificed dataset, generated using the corresponding script in this repo
+
+Author: Natasha Clarke; last edit 2024-07-16
+
+"""
+
+import argparse
+import pandas as pd
+from pathlib import Path
+
+
+def create_master_df(root_p, datasets):
+ # Create one df of QC results for specified datasets
+ cols = [
+ "identifier",
+ "participant_id",
+ "ses",
+ "run",
+ "age",
+ "sex",
+ "site",
+ "diagnosis",
+ "dataset",
+ "mean_fd_raw",
+ "mean_fd_scrubbed",
+ "proportion_kept",
+ "functional_dice",
+ "pass_func_qc",
+ "anatomical_dice",
+ "pass_anat_qc",
+ "pass_all_qc",
+ "different_func_affine",
+ ] # scanner
+
+ master_df = pd.DataFrame()
+ for dataset in datasets:
+ path_template = "wrangling-phenotype/outputs/{dataset}_qc_pheno.tsv"
+ df_p = root_p / path_template.format(dataset=dataset)
+ df = pd.read_csv(df_p, sep="\t", dtype={"participant_id": str})
+ df = df[cols]
+ master_df = pd.concat([master_df, df], ignore_index=True)
+
+ return master_df
+
+
+def filter_diagnoses(df, diagnoses):
+ return df[df["diagnosis"].isin(diagnoses)]
+
+
+def summarise_qc(df):
+ # Across all scans i.e. rows
+ scan_summary = (
+ master_df.groupby("dataset")
+ .agg(
+ total_scans=("participant_id", "count"), # Count all scans per dataset
+ total_scan_passes=(which_qc_col, "sum"), # Sum passes per dataset
+ )
+ .reset_index()
+ )
+
+ # Across sessions i.e. only count one run per session
+ session_passes = (
+ master_df.groupby(["dataset", "participant_id", "ses"])[which_qc_col]
+ .any()
+ .reset_index()
+ )
+ session_summary = session_passes.groupby("dataset").agg(
+ total_sessions=("ses", "count"), # Count unique sessions per dataset
+ total_session_passes=(which_qc_col, "sum"), # Count sessions that passed
+ )
+
+ # Across participants i.e. only count one session per participant
+ subject_summary = (
+ session_passes.groupby(["dataset", "participant_id"])[which_qc_col]
+ .any()
+ .reset_index()
+ )
+ subject_summary = subject_summary.groupby("dataset").agg(
+ total_subjects=("participant_id", "count"), # Count unique subjects per dataset
+ total_subject_passes=(which_qc_col, "sum"), # Count subjects that passed
+ )
+
+ return (scan_summary, session_summary, subject_summary)
+
+
+def combine_summaries(scan_summary, session_summary, subject_summary):
+ # Ensure 'dataset' is a column for joining
+ session_summary = session_summary.reset_index()
+ subject_summary = subject_summary.reset_index()
+
+ # Merge summaries
+ combined_summary = pd.merge(
+ scan_summary, session_summary, on="dataset", how="outer"
+ )
+ combined_summary = pd.merge(
+ combined_summary, subject_summary, on="dataset", how="outer"
+ )
+
+ combined_summary.columns = [
+ "dataset",
+ "Total Scans",
+ "Total Scan Passes",
+ "Total Sessions",
+ "Total Session Passes",
+ "Total Subjects",
+ "Total Subject Passes",
+ ]
+
+ return combined_summary
+
+
+if __name__ == "__main__":
+ parser = argparse.ArgumentParser(
+ description="Generate QC summary metrics and .tsv of scans passed QC"
+ )
+ parser.add_argument("--root_p", type=Path, help="Root path for data")
+ parser.add_argument(
+ "--datasets", nargs="+", type=str, help="List of datasets to process"
+ )
+ parser.add_argument(
+ "--which_qc_col",
+ type=str,
+ help="Column for evaluating QC. Default=pass_func_qc",
+ )
+ parser.add_argument(
+ "--diagnoses", nargs="+", type=str, help="Diagnoses of interest"
+ )
+
+ args = parser.parse_args()
+ datasets = args.datasets
+ which_qc_col = args.which_qc_col or "pass_func_qc"
+ diagnoses = args.diagnoses if args.diagnoses else None
+ output_p = args.root_p / "wrangling-phenotype/outputs"
+
+ # Create df of QC results and pheno data for specified datasets
+ master_df = create_master_df(args.root_p, datasets)
+
+ # Optionally filter df for specific diagnoses
+ if diagnoses:
+ master_df = filter_diagnoses(master_df, diagnoses)
+
+ # Summarise QC results
+ # Fill in blank session labels with a dummy variable to enable aggregation
+ master_df["ses"] = master_df["ses"].fillna("dummy").replace("", "dummy")
+ scan_summary, session_summary, subject_summary = summarise_qc(master_df)
+ qc_summary = combine_summaries(scan_summary, session_summary, subject_summary)
+ # Remove dummy variable
+ master_df["ses"].replace("dummy", "", inplace=True)
+
+ # Filter for scans that passed QC based on specified column (default = pass_func_qc)
+ filtered_df = master_df[master_df[which_qc_col] == True].copy()
+
+ # Save output
+ qc_summary.to_csv(output_p / "qc_summary.tsv", sep="\t", index=False)
+ filtered_df.to_csv(output_p / "passed_qc_master.tsv", sep="\t", index=False)
+
+ print(f"Data have been processed and output to {output_p}")
diff --git a/code/filtered_qc_report.ipynb b/code/filtered_qc_report.ipynb
new file mode 100644
index 0000000..2faf0ea
--- /dev/null
+++ b/code/filtered_qc_report.ipynb
@@ -0,0 +1,556 @@
+{
+ "cells": [
+ {
+ "cell_type": "code",
+ "execution_count": 1,
+ "id": "0a7e29ba",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "import matplotlib.pyplot as plt\n",
+ "import seaborn as sns\n",
+ "import pandas as pd\n",
+ "\n",
+ "from pathlib import Path"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 2,
+ "id": "2ddd34a9",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def boxplot_diagnosis_dataset(df, y_column):\n",
+ " plt.figure(figsize=(12, 8))\n",
+ " sns.boxplot(x='diagnosis', y=y_column, hue='dataset', data=df, palette=\"colorblind\")\n",
+ " plt.title(f'{y_column} by diagnosis across datasets')\n",
+ " plt.xlabel('Diagnosis')\n",
+ " plt.ylabel(y_column)\n",
+ " plt.xticks(rotation=45)\n",
+ " plt.legend(title='dataset')\n",
+ " plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 3,
+ "id": "26a1e863",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def histogram_diagnosis_dataset(df, y_column):\n",
+ " g = sns.FacetGrid(df, col=\"dataset\", row=\"diagnosis\", hue=\"diagnosis\", \n",
+ " height=4, sharex=True, sharey=False, palette=\"colorblind\")\n",
+ "\n",
+ " # Map the plt.hist function to each subplot\n",
+ " g.map(plt.hist, y_column, bins=10, edgecolor=\"k\").add_legend()\n",
+ "\n",
+ " g.fig.suptitle(f'{y_column} by diagnosis across datasets', fontsize=16, y=1.05)\n",
+ " g.set_axis_labels(y_column, 'Count')\n",
+ " g.set_titles(\"{col_name}\", size=14)\n",
+ " plt.show()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 4,
+ "id": "00946003",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "def boxplot_diagnosis_sex(df, y_column, group):\n",
+ " plt.figure(figsize=(6, 4))\n",
+ " sns.boxplot(data=df[df['diagnosis'] == group], x='sex', y=y_column, palette=\"colorblind\")\n",
+ " plt.title(f'{y_column} by diagnosis and sex ({group} only)')\n",
+ " plt.xlabel('Sex')\n",
+ " plt.ylabel(y_column)\n",
+ " plt.xticks(rotation=45)\n",
+ " plt.show()"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "7ce11f4d",
+ "metadata": {},
+ "source": [
+ "# Load TSV of scans passed QC, generated using script filter_pheno_qc.py"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 5,
+ "id": "981dd01c",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "passed_qc_p = Path(\"/home/neuromod/wrangling-phenotype/outputs/passed_qc_master.tsv\")\n",
+ "df = pd.read_csv(passed_qc_p, sep=\"\\t\")"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 6,
+ "id": "3ea537a0",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "dataset\n",
+ "adni 3294\n",
+ "cimaq 268\n",
+ "cobre 120\n",
+ "ds000030 154\n",
+ "hcpep 689\n",
+ "oasis3 6913\n",
+ "srpbs 711\n",
+ "dtype: int64"
+ ]
+ },
+ "execution_count": 6,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# How many scans passed QC?\n",
+ "df.groupby('dataset').size()"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 7,
+ "id": "211de586",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/plain": [
+ "dataset\n",
+ "adni 501\n",
+ "cimaq 107\n",
+ "cobre 120\n",
+ "ds000030 154\n",
+ "hcpep 177\n",
+ "oasis3 783\n",
+ "srpbs 711\n",
+ "Name: participant_id, dtype: int64"
+ ]
+ },
+ "execution_count": 7,
+ "metadata": {},
+ "output_type": "execute_result"
+ }
+ ],
+ "source": [
+ "# How many individual participants passed QC?\n",
+ "df.groupby('dataset').nunique()['participant_id']"
+ ]
+ },
+ {
+ "cell_type": "markdown",
+ "id": "28ab4a0a",
+ "metadata": {},
+ "source": [
+ "# Take a look at diagnoses"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 8,
+ "id": "695cfe74",
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "# Optionally re-classify some diagnoses, for ease\n",
+ "df[\"diagnosis\"] = df[\"diagnosis\"].replace({\"EMCI\": \"MCI\", \"LMCI\": \"MCI\", \"PSYC\":\"SCHZ\", \"ADD(M)\": \"ADD\"})"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": 9,
+ "id": "e0e99a72",
+ "metadata": {},
+ "outputs": [
+ {
+ "data": {
+ "text/html": [
+ "\n",
+ "\n",
+ "
\n",
+ " \n",
+ " \n",
+ " | dataset | \n",
+ " adni | \n",
+ " cimaq | \n",
+ " cobre | \n",
+ " ds000030 | \n",
+ " hcpep | \n",
+ " oasis3 | \n",
+ " srpbs | \n",
+ " TOTAL | \n",
+ "
\n",
+ " \n",
+ " | diagnosis | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ " | \n",
+ "
\n",
+ " \n",
+ " \n",
+ " \n",
+ " | ADD | \n",
+ " 68 | \n",
+ " 3 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 156 | \n",
+ " 0 | \n",
+ " 227 | \n",
+ "
\n",
+ " \n",
+ " | CON | \n",
+ " 309 | \n",
+ " 59 | \n",
+ " 67 | \n",
+ " 110 | \n",
+ " 57 | \n",
+ " 647 | \n",
+ " 605 | \n",
+ " 1837 | \n",
+ "
\n",
+ " \n",
+ " | MCI | \n",
+ " 163 | \n",
+ " 51 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 4 | \n",
+ " 0 | \n",
+ " 218 | \n",
+ "
\n",
+ " \n",
+ " | SCHZ | \n",
+ " 0 | \n",
+ " 0 | \n",
+ " 53 | \n",
+ " 44 | \n",
+ " 120 | \n",
+ " 0 | \n",
+ " 106 | \n",
+ " 323 | \n",
+ "
\n",
+ " \n",
+ " | TOTAL | \n",
+ " 501 | \n",
+ " 107 | \n",
+ " 120 | \n",
+ " 154 | \n",
+ " 177 | \n",
+ " 783 | \n",
+ " 711 | \n",
+ " 2498 | \n",
+ "
\n",
+ " \n",
+ "
\n",
+ "