Fix up code.

Author: tsalo

analysis/compare_tedana_aroma.py

@@ -0,0 +1,142 @@
+"""Compare the results of TEDANA and AROMA."""
+
+import os
+from glob import glob
+
+import matplotlib.pyplot as plt
+import pandas as pd
+import seaborn as sns
+
+
+if __name__ == "__main__":
+    in_dir = "/cbica/projects/pafin/derivatives/tedana+aroma"
+    out_dir = "../figures"
+
+    files = sorted(
+        glob(
+            os.path.join(
+                in_dir,
+                "sub-*",
+                "ses-1",
+                "func",
+                "*_desc-tedana+aroma_metrics.tsv",
+            )
+        )
+    )
+    group_df = pd.DataFrame(
+        columns=[
+            "file",
+            "subject",
+            "denoising",
+            "n_components",
+            "n_accepted",
+            "n_rejected_aroma",
+            "n_rejected_tedana",
+            "n_rejected_both",
+            "varex_accepted",
+            "varex_rejected_aroma",
+            "varex_rejected_tedana",
+            "varex_rejected_both",
+            "varex_unmodeled",
+        ],
+    )
+    for file in files:
+        df = pd.read_table(file)
+        subject = os.path.basename(file).split("_")[0]
+        file_dict = {
+            "file": os.path.basename(file),
+            "subject": subject,
+            "denoising": "NORDIC" if "nordic" in os.path.basename(file) else "None",
+            "n_components": df.shape[0],
+            "n_accepted": df[df["classification"] == "accepted"].shape[0],
+            "varex_accepted": df[df["classification"] == "accepted"][
+                "total_variance_explained"
+            ].sum(),
+            "varex_unmodeled": 1 - df["total_variance_explained"].sum(),
+        }
+
+        rejected_df = df.loc[df["classification"] == "rejected"]
+        rejected_both_df = rejected_df.loc[
+            rejected_df["classification_tags"].str.contains("TEDANA")
+            & rejected_df["classification_tags"].str.contains("AROMA")
+        ]
+        rejected_aroma_df = rejected_df.loc[
+            rejected_df["classification_tags"].str.contains("AROMA")
+            & ~rejected_df["classification_tags"].str.contains("TEDANA")
+        ]
+        rejected_tedana_df = rejected_df.loc[
+            ~rejected_df["classification_tags"].str.contains("AROMA")
+            & rejected_df["classification_tags"].str.contains("TEDANA")
+        ]
+
+        file_dict["n_rejected_both"] = rejected_both_df.shape[0]
+        file_dict["n_rejected_aroma"] = rejected_aroma_df.shape[0]
+        file_dict["n_rejected_tedana"] = rejected_tedana_df.shape[0]
+        file_dict["varex_rejected_both"] = rejected_both_df["total_variance_explained"].sum()
+        file_dict["varex_rejected_aroma"] = rejected_aroma_df["total_variance_explained"].sum()
+        file_dict["varex_rejected_tedana"] = rejected_tedana_df["total_variance_explained"].sum()
+
+        group_df = group_df.append(file_dict, ignore_index=True)
+
+    group_df.to_csv("../data/bold_denoising_metrics.tsv", index=False, sep="\t")
+
+    # Boxplot of variance explained, organized as "accepted", "rejected by AROMA",
+    # "rejected by TEDANA", "rejected by both", "unmodeled" across runs
+    df_varex = group_df.melt(
+        id_vars=["subject", "denoising"],
+        value_vars=[
+            "varex_accepted",
+            "varex_rejected_aroma",
+            "varex_rejected_tedana",
+            "varex_rejected_both",
+            "varex_unmodeled",
+        ],
+        value_name="Variance Explained",
+        var_name="Classification",
+    )
+
+    sns.set_theme(style="ticks")
+    f, ax = plt.subplots(figsize=(7, 6))
+    sns.boxenplot(
+        df_varex,
+        x="Variance Explained",
+        y="Classification",
+        hue="denoising",
+        palette="vlag",
+    )
+    ax.xaxis.grid(True)
+    ax.set_xlim(0, 1)
+    ax.set(ylabel="")
+    sns.despine(trim=True, left=True)
+    f.savefig(os.path.join(out_dir, "bold_denoising_variance_explained.png"))
+    plt.close()
+
+    # Boxplot of number of components, organized as "accepted", "rejected by AROMA",
+    # "rejected by TEDANA", "rejected by both" across runs
+    df_ncomps = group_df.melt(
+        id_vars=["subject", "denoising"],
+        value_vars=[
+            "n_accepted",
+            "n_rejected_aroma",
+            "n_rejected_tedana",
+            "n_rejected_both",
+        ],
+        value_name="Number of Components",
+        var_name="Classification",
+    )
+
+    sns.set_theme(style="ticks")
+    f, ax = plt.subplots(figsize=(7, 6))
+    sns.boxenplot(
+        df_ncomps,
+        x="Number of Components",
+        y="Classification",
+        hue="denoising",
+        palette="vlag",
+    )
+    ax.xaxis.grid(True)
+    ax.set_xlim(0, None)
+    ax.set(ylabel="")
+    sns.despine(trim=True, left=True)
+    f.savefig(os.path.join(out_dir, "bold_denoising_n_components.png"))
+    plt.close()

processing/run_tedana.py

@@ -262,22 +262,35 @@ def run_tedana_aroma(raw_dir, aroma_dir, tedana_out_dir, tedana_aroma_out_dir):
         for i_row, aroma_row in aroma_df.iterrows():
             aroma_clf = aroma_row["classification"]
             tedana_rationale = tedana_df.loc[i_row, "classification_tags"]
+            tedana_rationales = tedana_rationale.split(";")
 
+            aroma_rationales = []
             if aroma_clf == "rejected":
                 aroma_rationales = aroma_row["rationale"].split(";")
-                aroma_rationales = [
-                    f"AROMA {aroma_rationale}" for aroma_rationale in aroma_rationales
-                ]
-                aroma_rationale = ";".join(aroma_rationales)
-                tedana_rationale += f";{aroma_rationale}"
-                rationales = tedana_rationale.split(";")
-                rationales = [
-                    r for r in rationales if r not in ["low variance", "accept borderline"]
+                aroma_rationales = [f"AROMA {rationale}" for rationale in aroma_rationales]
+
+                # Remove tags equivalent to "ignored" classification, for MIR
+                tedana_rationales = [
+                    r for r in tedana_rationales if r not in ["low variance", "accept borderline"]
                 ]
-                tedana_rationale = ";".join(rationales)
 
                 tedana_df.iloc[i_row, "classification"] = "rejected"
-                tedana_df.iloc[i_row, "classification_tags"] = tedana_rationale
+
+            tedana_rationales = [f"TEDANA {rationale}" for rationale in tedana_rationales]
+            rationales = tedana_rationales + aroma_rationales
+            tedana_df.iloc[i_row, "classification_tags"] = ";".join(rationales)
+
+            # Add other columns from aroma_df to tedana_df
+            other_cols = [
+                "edge_fract",
+                "csf_fract",
+                "max_RP_corr",
+                "HFC",
+                "model_variance_explained",
+                "total_variance_explained",
+            ]
+            for col in other_cols:
+                tedana_df.loc[i_row, col] = aroma_row[col]
 
         # Save the combined classifications
         tedana_aroma_run_out_dir = os.path.join(