Repository: ME-ICA/tedana

Issue:
Add global signal (if `--gscontrol gsr` is used) to component plots
### Summary
I think it might be useful to include the global signal in the HTML report, and I was thinking that it could be treated as a component for the sake of the report. It would also be nice to see how much variance the global signal explains.

Discussion:
@tsalo: Alternatively, we could just add brain plots and time series plots to the HTML reports in new sections, like we've recently done with the T2*/S0 maps and the adaptive mask.

We wouldn't get as good a sense of the amount of variance explained in that case though.@tsalo: I can plot the global signal map, but I wonder if it would make more sense to plot a variance explained map instead?

Solution patch:
diff --git a/docs/outputs.rst b/docs/outputs.rst
index 2ed3ddb..3c59412 100644
--- a/docs/outputs.rst
+++ b/docs/outputs.rst
@@ -618,6 +618,18 @@ relatively high RMSE values for runs or timepoints might be a marker of an under
   :height: 400px
 
 
+***************************
+Global Signal Removal Plots
+***************************
+
+Below the decay model fit plots are the global signal removal plots, if ``--gscontrol`` was used.
+These plots show the voxel-wise global signal weights and the global signal time series.
+
+.. image:: /_static/gsr_plots.png
+  :align: center
+  :height: 400px
+
+
 **************************
 Citable workflow summaries
 **************************
diff --git a/tedana/gscontrol.py b/tedana/gscontrol.py
index acd35cb..f17e7b0 100644
--- a/tedana/gscontrol.py
+++ b/tedana/gscontrol.py
@@ -120,7 +120,13 @@ def gscontrol_raw(
     # Project global signal (but not Legendre bases) out of optimally combined data
     betas = np.linalg.lstsq(glbase, data_optcom_masked.T, rcond=None)[0]
     gs_fitted = np.dot(glbase[:, :1], betas[:1, :]).T
-    data_optcom_nogs = data_optcom_masked - gs_fitted + temporal_mean
+    data_optcom_nogs = data_optcom_masked - gs_fitted
+
+    # Calculate the variance explained by the global signal
+    varexpl = 100 * (1 - ((data_optcom_nogs**2).sum() / (data_optcom_masked**2).sum()))
+    LGR.info(f"Variance in optimally combined data explained by global signal: {varexpl:.02f}%")
+
+    data_optcom_nogs += temporal_mean
     data_optcom_nogs = utils.unmask(data_optcom_nogs, temporal_mean_mask)
     io_generator.save_file(data_optcom_nogs, "removed gs combined img")
 
@@ -257,12 +263,18 @@ def minimum_image_regression(
     glsig_df = pd.DataFrame(data=gs_ts.T, columns=["mir_global_signal"])
     io_generator.add_df_to_file(glsig_df, "confounds tsv")
 
-    # Remove T1-like global signal from MEHK time series
-    mehk_no_t1_gs = mehk_ts - np.dot(
-        np.linalg.lstsq(gs_ts.T, mehk_ts.T, rcond=None)[0].T,
-        gs_ts,
+    # Calculate the T1-like global signal
+    t1_gs_data = np.dot(np.linalg.lstsq(gs_ts.T, mehk_ts.T, rcond=None)[0].T, gs_ts)
+
+    # Calculate the variance explained by the T1-like global signal
+    varexpl = 100 * (1 - (((data_optcom_z - t1_gs_data) ** 2).sum() / (data_optcom_z**2).sum()))
+    LGR.info(
+        f"Variance in optimally combined data explained by T1-like global signal: {varexpl:.02f}%"
     )
 
+    # Remove T1-like global signal from MEHK time series
+    mehk_no_t1_gs = mehk_ts - t1_gs_data
+
     # Make denoised version of T1-corrected time series
     medn_ts = optcom_mean + ((mehk_no_t1_gs + resid) * optcom_std)
     io_generator.save_file(utils.unmask(medn_ts, mask), "mir denoised img")
diff --git a/tedana/reporting/data/html/report_body_template.html b/tedana/reporting/data/html/report_body_template.html
index 75304e2..6c87888 100644
--- a/tedana/reporting/data/html/report_body_template.html
+++ b/tedana/reporting/data/html/report_body_template.html
@@ -252,6 +252,40 @@
   </div>
   {% endif %}
 </div>
+
+{% if gsr or mir %}
+<div class="carpet-plots">
+  <h2>Global Signal Control</h2>
+    {% if gsr %}
+      <h3>Global Signal Regression Weights</h3>
+      <div class="carpet-plots-image">
+        <img
+          id="gsrBrainPlot"
+          src="{{ gsrBrainPlot }}"
+          style="width: 95%; max-width: 1200px; height: auto"
+        />
+      </div>
+    {% endif %}
+    {% if mir %}
+      <h3>Minimum Image Regression Weights</h3>
+      <div class="carpet-plots-image">
+        <img
+          id="mirBrainPlot"
+          src="{{ mirBrainPlot }}"
+          style="width: 95%; max-width: 1200px; height: auto"
+        />
+      </div>
+    {% endif %}
+    <h3>Global Signal Timeseries</h3>
+    <div class="carpet-plots-image">
+      <img
+        id="gsrTimeseries"
+        src="{{ gsrTimeseries }}"
+        style="width: 95%; max-width: 1200px; height: auto"
+      />
+    </div>
+</div>
+{% endif %}
 <div class="info">
   <h1>Info</h1>
   {{ info }}
diff --git a/tedana/reporting/html_report.py b/tedana/reporting/html_report.py
index 4c50cfa..0db1590 100644
--- a/tedana/reporting/html_report.py
+++ b/tedana/reporting/html_report.py
@@ -183,6 +183,16 @@ def _update_template_bokeh(
     rmse_brain = f"./figures/{rmse_brain_filename}"
     rmse_timeseries = f"./figures/{rmse_timeseries_filename}"
 
+    # Check for gscontrol images
+    gsr_brain_filename = f"{prefix}gsr_boldmap.svg"
+    mir_brain_filename = f"{prefix}T1likeEffect_boldmap.svg"
+    gsr_timeseries_filename = f"{prefix}gscontrol_bold.svg"
+    gsr_brain = f"./figures/{gsr_brain_filename}"
+    mir_brain = f"./figures/{mir_brain_filename}"
+    gsr_timeseries = f"./figures/{gsr_timeseries_filename}"
+    gsr = gsr_brain_filename in files_in_figures and gsr_timeseries_filename in files_in_figures
+    mir = mir_brain_filename in files_in_figures and gsr_timeseries_filename in files_in_figures
+
     # Check if each set of images exists
     t2star_exists = (
         t2star_brain_filename in files_in_figures and t2star_histogram_filename in files_in_figures
@@ -224,6 +234,11 @@ def _update_template_bokeh(
         rmseBrainPlot=rmse_brain,
         rmseTimeseries=rmse_timeseries,
         rmseExists=rmse_exists,
+        gsr=gsr,
+        mir=mir,
+        gsrBrainPlot=gsr_brain,
+        mirBrainPlot=mir_brain,
+        gsrTimeseries=gsr_timeseries,
         references=references,
         javascript=bokeh_js,
         buttons=buttons,
diff --git a/tedana/reporting/static_figures.py b/tedana/reporting/static_figures.py
index 6f8682b..16f6273 100644
--- a/tedana/reporting/static_figures.py
+++ b/tedana/reporting/static_figures.py
@@ -8,6 +8,7 @@ from io import BytesIO
 import matplotlib
 import nibabel as nb
 import numpy as np
+import scipy.stats as sstats
 
 matplotlib.use("AGG")
 import matplotlib.pyplot as plt
@@ -814,3 +815,92 @@ def plot_adaptive_mask(
     )
     adaptive_mask_plot = f"{io_generator.prefix}adaptive_mask.svg"
     fig.savefig(os.path.join(io_generator.out_dir, "figures", adaptive_mask_plot))
+
+
+def plot_gscontrol(
+    *,
+    io_generator: io.OutputGenerator,
+    gscontrol: list,
+):
+    """Plot the results of the gscontrol steps.
+
+    Parameters
+    ----------
+    io_generator : :obj:`~tedana.io.OutputGenerator`
+        The output generator for this workflow.
+    gscontrol : list
+        List of gscontrol methods applied.
+    """
+    import pandas as pd
+
+    if "gsr" in gscontrol:
+        gsr_img = io_generator.get_name("gs img")
+        gsr_plot = f"{io_generator.prefix}gsr_boldmap.svg"
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                "ignore",
+                message="A non-diagonal affine.*",
+                category=UserWarning,
+            )
+            plotting.plot_stat_map(
+                gsr_img,
+                bg_img=None,
+                display_mode="mosaic",
+                symmetric_cbar=True,
+                black_bg=True,
+                cmap="coolwarm",
+                annotate=False,
+                output_file=os.path.join(io_generator.out_dir, "figures", gsr_plot),
+            )
+
+    if "mir" in gscontrol:
+        mir_img = io_generator.get_name("t1 like img")
+        mir_plot = f"{io_generator.prefix}T1likeEffect_boldmap.svg"
+        with warnings.catch_warnings():
+            warnings.filterwarnings(
+                "ignore",
+                message="A non-diagonal affine.*",
+                category=UserWarning,
+            )
+            plotting.plot_stat_map(
+                mir_img,
+                bg_img=None,
+                display_mode="mosaic",
+                symmetric_cbar=True,
+                black_bg=True,
+                cmap="coolwarm",
+                annotate=False,
+                output_file=os.path.join(io_generator.out_dir, "figures", mir_plot),
+            )
+
+    if "gsr" in gscontrol or "mir" in gscontrol:
+        confounds_file = io_generator.get_name("confounds tsv")
+        confounds_df = pd.read_table(confounds_file)
+
+        # Get repetition time from reference image
+        tr = io_generator.reference_img.header.get_zooms()[-1]
+        time_arr = np.arange(confounds_df.shape[0]) * tr
+
+        fig, ax = plt.subplots(figsize=(10, 6))
+        if "gsr" in gscontrol:
+            gs = confounds_df["global_signal"].values
+            gs_z = sstats.zscore(gs)
+            ax.plot(time_arr, gs_z, label="GSR", color="red")
+
+        if "mir" in gscontrol:
+            mir = confounds_df["mir_global_signal"].values
+            mir_z = sstats.zscore(mir)
+            ax.plot(time_arr, mir_z, label="MIR", color="blue")
+
+        ax.legend()
+        ax.set_title("Global Signals")
+        ax.set_xlabel("Time (s)")
+        ax.set_ylabel("Signal (z-scored)")
+        ax.set_xlim(0, time_arr[-1])
+        fig.savefig(
+            os.path.join(
+                io_generator.out_dir,
+                "figures",
+                f"{io_generator.prefix}gscontrol_bold.svg",
+            )
+        )
diff --git a/tedana/workflows/tedana.py b/tedana/workflows/tedana.py
index 1406083..bea2b57 100644
--- a/tedana/workflows/tedana.py
+++ b/tedana/workflows/tedana.py
@@ -1136,6 +1136,12 @@ def tedana_workflow(
                 adaptive_mask=masksum_denoise,
             )
 
+        if gscontrol:
+            reporting.static_figures.plot_gscontrol(
+                io_generator=io_generator,
+                gscontrol=gscontrol,
+            )
+
         LGR.info("Generating dynamic report")
         reporting.generate_report(io_generator, cluster_labels, similarity_t_sne)
 


Test patch:
diff --git a/tedana/tests/data/fiu_four_echo_outputs.txt b/tedana/tests/data/fiu_four_echo_outputs.txt
index 9cb0832..bdf1639 100644
--- a/tedana/tests/data/fiu_four_echo_outputs.txt
+++ b/tedana/tests/data/fiu_four_echo_outputs.txt
@@ -70,6 +70,9 @@ figures/sub-01_carpet_rejected.svg
 figures/sub-01_carpet_optcom_nogsr.svg
 figures/sub-01_carpet_denoised_mir.svg
 figures/sub-01_carpet_accepted_mir.svg
+figures/sub-01_T1likeEffect_boldmap.svg
+figures/sub-01_gscontrol_bold.svg
+figures/sub-01_gsr_boldmap.svg
 figures/sub-01_comp_000.png
 figures/sub-01_comp_001.png
 figures/sub-01_comp_002.png