From 46b699d9875f55d0a7cebdaad06bf94f5b2fd96a Mon Sep 17 00:00:00 2001
From: Taylor Salo <salot@pennmedicine.upenn.edu>
Date: Tue, 19 Nov 2024 12:30:40 -0500
Subject: [PATCH] Make connectivity field names compatible with MATLAB (#166)

* Add "atlas_" prefix to connectivity file fields.

* Sanitize measure names.

* Update connectivity_matrices.rst

* Update connectivity_matrices.rst
---
 docs/connectivity_matrices.rst    | 40 ++++++++++++++++++++++++++++++-
 qsirecon/interfaces/dsi_studio.py | 10 ++++----
 qsirecon/interfaces/mrtrix.py     |  6 ++---
 3 files changed, 48 insertions(+), 8 deletions(-)

diff --git a/docs/connectivity_matrices.rst b/docs/connectivity_matrices.rst
index c3c75bab..7a21d386 100644
--- a/docs/connectivity_matrices.rst
+++ b/docs/connectivity_matrices.rst
@@ -14,7 +14,7 @@ these pipelines and only look at those numbers.
 If you look at more than one weighting method be sure to adjust your statistics for the
 additional comparisons.
 
-To skip this step in your workflow, you can modify an existing recon pipeline by removing the 
+To skip this step in your workflow, you can modify an existing recon pipeline by removing the
 ``action: connectivity`` section from the yaml file.
 
 .. _connectivity_atlases:
@@ -97,6 +97,44 @@ It is essential that your images are in the LPS+ orientation and have the sform
 **Be sure to check for alignment and orientation** in your outputs.
 
 
+*********************
+Connectivity Measures
+*********************
+
+Connectivity measures are bundled together in binary ``.mat`` files,
+rather than as atlas- and measure-specific tabular files.
+
+.. warning::
+
+   We ultimately plan to organize the connectivity matrices accoring to the BIDS-Connectivity BEP,
+   wherein each measure from each atlas is stored in a separate file.
+
+   Therefore, this organization will change in the future.
+
+.. code-block::
+
+   qsirecon/
+      derivatives/
+         qsirecon-<suffix>/
+            sub-<label>/[ses-<label>/]
+               dwi/
+                  <source_entities>_connectivity.mat
+
+The ``.mat`` file contains a dictionary with all of the connectivity measures specified
+by the recon spec for all of the different atlases specified by the user.
+
+For example, in the case where a user has selected a single atlas (``<atlas>``) and
+the recon spec specifies a single connectivity measure (``<measure>``),
+the ``.mat`` file will contain the following keys:
+
+.. code-block::
+
+   command                                # The command that was run
+   atlas_<atlas>_region_ids               # The region ids for the atlas (1 x n_parcels array)
+   atlas_<atlas>_region_labels            # The region labels for the atlas (1 x n_parcels array)
+   atlas_<atlas>_<measure>_connectivity   # The connectivity matrix for the atlas and measure (n_parcels x n_parcels array)
+
+
 **********
 References
 **********
diff --git a/qsirecon/interfaces/dsi_studio.py b/qsirecon/interfaces/dsi_studio.py
index 7010f38f..2b3a7520 100644
--- a/qsirecon/interfaces/dsi_studio.py
+++ b/qsirecon/interfaces/dsi_studio.py
@@ -334,15 +334,15 @@ def _post_run_hook(self, runtime):
         # Aggregate the connectivity/network data from DSI Studio
         official_labels = atlas_labels_df["index"].values.astype(int)
         connectivity_data = {
-            f"{atlas_name}_region_ids": official_labels,
-            f"{atlas_name}_region_labels": atlas_labels_df["label"].values,
+            f"atlas_{atlas_name}_region_ids": official_labels,
+            f"atlas_{atlas_name}_region_labels": atlas_labels_df["label"].values,
         }
 
         # Gather the connectivity matrices
         matfiles = glob(op.join(runtime.cwd, "*.connectivity.mat"))
         for matfile in matfiles:
             measure = "_".join(matfile.split(".")[-4:-2])
-            connectivity_data[f"{atlas_name}_{measure}_connectivity"] = (
+            connectivity_data[f"atlas_{atlas_name}_{measure}_connectivity"] = (
                 _sanitized_connectivity_matrix(matfile, official_labels)
             )
 
@@ -552,7 +552,9 @@ def _sanitized_network_measures(network_txt, official_labels, atlas_name, measur
     assert np.all(truncated_labels == network_region_ids)
 
     for net_measure_name, net_measure_data in network_data.items():
-        variable_name = atlas_name + "_" + measure + "_" + net_measure_name
+        net_measure_name = net_measure_name.replace("-", "_")
+        measure_name = measure.replace("-", "_")
+        variable_name = f"atlas_{atlas_name}_{measure_name}_{net_measure_name}"
         if type(net_measure_data) is np.ndarray:
             tmp = np.zeros(n_atlas_labels)
             tmp[in_this_mask] = net_measure_data
diff --git a/qsirecon/interfaces/mrtrix.py b/qsirecon/interfaces/mrtrix.py
index faf860ef..d9cabd78 100644
--- a/qsirecon/interfaces/mrtrix.py
+++ b/qsirecon/interfaces/mrtrix.py
@@ -744,13 +744,13 @@ def _post_run_hook(self, runtime):
 
         # Aggregate the connectivity/network data from DSI Studio
         connectivity_data = {
-            f"{atlas_name}_region_ids": atlas_labels_df["index"].values.astype(int),
-            f"{atlas_name}_region_labels": atlas_labels_df["label"].values,
+            f"atlas_{atlas_name}_region_ids": atlas_labels_df["index"].values.astype(int),
+            f"atlas_{atlas_name}_region_labels": atlas_labels_df["label"].values,
         }
 
         # get the connectivity matrix
         prefix = f"{atlas_name}_{self.inputs.measure}"
-        connectivity_data[f"{prefix}_connectivity"] = np.loadtxt(
+        connectivity_data[f"atlas_{prefix}_connectivity"] = np.loadtxt(
             self.inputs.out_file, delimiter=","
         )
         connectivity_data["command"] = self.cmdline