High-level Acetylation Comparisons

aml_outcome_disparities/analysis/fig3/figure3a.py

@@ -0,0 +1,92 @@
+"""Plots figure 3a: S100A8 Acetylation vs. Race."""
+
+from matplotlib.patches import Patch
+import numpy as np
+import pandas as pd
+
+from pilot.data_import import import_acetyl, import_meta, syn_login
+from pilot.figures.figure_setup import compare_means, get_setup
+
+
+def make_figure():
+    # Import meta-data, acetyl
+    syn = syn_login()
+    meta = import_meta(syn)
+    acetyl = pd.concat(import_acetyl(syn))
+
+    # Trim to S100A8 sites
+    acetyl = acetyl.loc[
+        :,
+        acetyl.columns.str.startswith("S100A8")
+    ]
+
+    # Trim to only Black and White patients
+    race = meta.loc[acetyl.index, "Race"]
+    acetyl = acetyl.loc[race.isin(["White", "Black"]), :]
+    race = race.loc[acetyl.index]
+
+    # Setup figure
+    fig, ax = get_setup(1, 1, fig_params={"figsize": (2, 3)})
+
+    # Iterate through each acetyl site, compare Black and White patients
+    for x_index, acetyl_site in enumerate(acetyl.columns):
+        acetyl_col = acetyl.loc[:, acetyl_site].dropna()
+        _race = race.loc[acetyl_col.index]
+        white_bp = ax.boxplot(
+            acetyl_col.loc[_race == "White"],
+            positions=[3 * x_index],
+            widths=0.8,
+            patch_artist=True,
+            notch=True,
+            medianprops={"color": "black"},
+            flierprops={"color": "tab:purple"},
+        )
+        black_bp = ax.boxplot(
+            acetyl_col.loc[_race == "Black"],
+            positions=[3 * x_index + 1],
+            widths=0.8,
+            patch_artist=True,
+            notch=True,
+            medianprops={"color": "black"},
+            flierprops={"color": "tab:red"},
+        )
+        white_bp["boxes"][0].set_facecolor("tab:purple")
+        black_bp["boxes"][0].set_facecolor("tab:red")
+
+    # Compare means between acetyl sites
+    ax = compare_means(
+        acetyl,
+        race,
+        1,
+        3,
+        ax,
+        star_only=True,
+        alternative="greater"
+    )
+
+    # Reformat ticks and labels
+    ax.set(
+        xticks=np.arange(0.5, 3 * acetyl.shape[1], 3),
+        ylabel="Log2 Expression",
+        ylim=(-8, 8)
+    )
+    ax.set_xticklabels(
+        acetyl.columns, rotation=45, ha="right", ma="right", va="top"
+    )
+
+    # Add legend for Black & White patients
+    legend_elements = [
+        Patch(
+            facecolor="tab:purple",
+            edgecolor="tab:purple",
+            label="White Patients",
+        ),
+        Patch(facecolor="tab:red", edgecolor="tab:red", label="Black Patients")
+    ]
+    ax.legend(handles=legend_elements, loc="lower left")
+
+    return fig
+
+
+if __name__ == "__main__":
+    make_figure()

aml_outcome_disparities/analysis/fig3/figure3b.py

@@ -0,0 +1,85 @@
+"""Plots figure 3b: S100A8 Acetylation & ILKAP."""
+
+from decimal import Decimal
+
+import matplotlib.pyplot as plt
+import numpy as np
+import pandas as pd
+
+from pilot.data_import import (import_acetyl, import_global, import_rna,
+                               syn_login)
+from pilot.figures.figure_setup import get_setup, run_ols
+
+GENE = "ILKAP"
+
+
+def make_figure():
+    syn = syn_login()
+    acetyl = pd.concat(import_acetyl(syn, global_correct=False))
+    prot = pd.concat(import_global(syn))
+    rna = np.log2(pd.concat(import_rna(syn)) + 1)
+
+    prot = prot.loc[
+        prot.index.intersection(acetyl.index),
+        :
+    ]
+    prot = prot.loc[
+        prot.index.intersection(rna.index),
+        :
+    ]
+    rna = rna.loc[prot.index, :]
+    acetyl = acetyl.loc[
+        prot.index,
+        acetyl.columns.str.startswith("S100A8")
+    ]
+
+    fig, axes = get_setup(
+        1,
+        2,
+        fig_params={
+            "figsize": (3.5, 3),
+            "width_ratios": [3.375, 0.125]
+        }
+    )
+
+    ax = axes[0]
+    colorbar_ax = axes[1]
+
+    scatter = ax.scatter(
+        rna.loc[:, GENE],
+        prot.loc[:, GENE],
+        s=3,
+        c=acetyl.loc[
+            :,
+            acetyl.columns.str.startswith("S100A8")
+        ].mean(axis=1),
+        cmap="coolwarm"
+    )
+    ols = run_ols(rna.loc[:, GENE], prot.loc[:, GENE], ax)
+
+    ax.text(
+        0.99,
+        0.01,
+        ha="right",
+        ma="right",
+        va="bottom",
+        s=f"R-squared: {round(ols.rsquared_adj, 3)}\nCoefficient p-value: "
+          f"{'{:.2E}'.format(Decimal(ols.pvalues.iloc[0]))}",
+        transform=ax.transAxes
+    )
+
+    plt.colorbar(
+        scatter,
+        label="S100A8 Acetylation",
+        cax=colorbar_ax
+    )
+    ax.set(
+        xlabel=f"RNA: {GENE}",
+        ylabel=f"Protein: {GENE}"
+    )
+
+    return fig
+
+
+if __name__ == "__main__":
+    make_figure()

aml_outcome_disparities/analysis/fig3/figure3c.py

@@ -0,0 +1,78 @@
+"""Plots figure 3c: S100A8 Acetylation vs. ILKAP Gap."""
+
+from decimal import Decimal
+
+import numpy as np
+import pandas as pd
+
+from pilot.data_import import (import_acetyl, import_global, import_rna,
+                               syn_login)
+from pilot.figures.figure_setup import get_setup, run_ols
+
+GENE = "ILKAP"
+
+
+def make_figure():
+    # Import meta-data, data, log-transform TPM
+    syn = syn_login()
+    acetyl = pd.concat(import_acetyl(syn, global_correct=False))
+    prot = pd.concat(import_global(syn))
+    rna = np.log2(pd.concat(import_rna(syn)) + 1)
+
+    # Trim to patients with all measurements
+    prot = prot.loc[
+        prot.index.intersection(acetyl.index),
+        :
+    ]
+    prot = prot.loc[prot.index.intersection(rna.index), :]
+    rna = rna.loc[prot.index, GENE]
+    diff = prot.loc[:, GENE] - rna
+
+    # Trim to S100A8 sites
+    acetyl = acetyl.loc[
+        prot.index,
+        acetyl.columns.str.startswith("S100A8")
+    ].mean(axis=1).dropna()
+
+    # Setup figure
+    fig, ax = get_setup(
+        1,
+        1,
+        fig_params={
+            "figsize": (3, 3)
+        }
+    )
+
+    # Plot acetylation vs. global/RNA differences
+    ax.scatter(
+        acetyl,
+        diff.loc[acetyl.index],
+        s=3
+    )
+
+    # Add best fit line
+    ols = run_ols(acetyl, diff.loc[acetyl.index], ax)
+
+    # Write OLS results
+    ax.text(
+        0.99,
+        0.01,
+        ha="right",
+        ma="right",
+        va="bottom",
+        s=f"R-squared: {round(ols.rsquared_adj, 3)}\nCoefficient p-value: "
+          f"{'{:.2E}'.format(Decimal(ols.pvalues.iloc[0]))}",
+        transform=ax.transAxes
+    )
+
+    # Label axes
+    ax.set(
+        xlabel=f"S100A8 Acetylation",
+        ylabel=f"{GENE}: Global - RNA"
+    )
+
+    return fig
+
+
+if __name__ == "__main__":
+    make_figure()