Update plots for reporting

Author: RJMW

beams/__main__.py

@@ -203,9 +203,12 @@ def main():
     parser_sr.add_argument('-i', '--intensity-matrix', type=str, required=False,
                            help="Tab-delimited intensity matrix.")
 
-    parser_sr.add_argument('-o', '--output', type=str, required=False,
+    parser_sr.add_argument('-o', '--output', type=str, required=True,
                            help="Output file for the summary")
 
+    parser_sr.add_argument('-p', '--pdf', type=str, required=False,
+                           help="Output pdf file for the summary plots")
+
     parser_sr.add_argument('-d', '--db', type=str, required=True,
                            help="Sqlite database that contains the results from the previous steps.")
 
@@ -298,6 +301,10 @@ def main():
         df = in_out.combine_peaklist_matrix(args.peaklist, args.intensity_matrix)
         df_out = annotation.summary(df, db=args.db, single_row=args.single_row, single_column=args.single_column, convert_rt=args.convert_rt, ndigits_mz=args.ndigits_mz)
         df_out.to_csv(args.output, sep=separators[args.sep], index=False)
+        if args.pdf:
+            plots.report(db=args.db, pdf_out=args.pdf,
+                         column_corr="r_value", column_pvalue="pvalue",
+                         column_ppm_error="ppm_error", column_adducts="adduct")
 
     if args.step == "start-gui":
         from PyQt5 import QtWidgets

beams/gui.py

@@ -490,8 +490,9 @@ def run(self):
 
             separators = {"tab": "\t", "comma": ","}
             df_out.to_csv(self.lineEdit_summary_filename.text(), sep=separators[self.comboBox_separator.currentText()], index=False)
-            plots.report(df=df_out, fn_pdf=os.path.join(os.path.dirname(self.lineEdit_summary_filename.text()), "report.pdf"),
-                         column_ppm_error="ppm_error", column_adducts="adduct")
+            ext = os.path.splitext(self.lineEdit_summary_filename.text())[1]
+            plots.report(db=self.lineEdit_sql_database.text(), pdf_out=str(self.lineEdit_summary_filename.text()).replace(ext, ".pdf"),
+                         column_corr="r_value", column_pvalue="p_value", column_ppm_error="ppm_error", column_adducts="adduct")
             print("Done")
             print("")
 

beams/plots.py

@@ -25,22 +25,102 @@
 import seaborn as sns
 
 
-def report(df, column_ppm_error, column_adducts, fn_pdf):
+import sys
+import sqlite3
+import numpy as np
+import pandas as pd
+import matplotlib
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+import seaborn as sns
+from matplotlib.backends.backend_pdf import PdfPages
+
+
+import matplotlib.pyplot as plt
+import matplotlib.gridspec as gridspec
+import seaborn as sns
+from matplotlib.backends.backend_pdf import PdfPages
+
+
+def plot_correlations(column_corr, column_pvalue, df):
+
+    fig = plt.figure(figsize=(8, 8))
+    fig.set_size_inches(8.27, 11.69)
+
+    gs = gridspec.GridSpec(3, 3)
+    ax_main = plt.subplot(gs[1:3, :2])
+    ax_x_dist = plt.subplot(gs[0, :2], sharex=ax_main)
+    ax_y_dist = plt.subplot(gs[1:3, 2], sharey=ax_main)
+
+    ax_main.grid(linestyle='dashed')
+    ax_x_dist.grid(linestyle='dashed')
+    ax_y_dist.grid(linestyle='dashed')
+
+    ax_main.set_axisbelow(True)
+    ax_x_dist.set_axisbelow(True)
+    ax_y_dist.set_axisbelow(True)
+
+    hb = ax_main.hexbin(x=column_corr, y=column_pvalue, data=df, gridsize=(20, 20), mincnt=1, extent=[-1, 1.0, 0, 0.1])
+    ax_main.set(xlabel="Correlation coefficient (R)", ylabel="P-value", xticks=np.arange(-1, 1.1, 0.1), yticks=np.arange(0.0,  0.105, 0.005))
+
+    ax_x_dist.hist(x=column_corr, data=df, bins=40, align='mid', color="lightblue")
+    ax_x_dist.set(ylabel='count', xlim=(-1, 1))
+    ax_x_dist.axvline(0, color='k', linestyle='dashed', linewidth=1)
+
+    ax_xcum_dist = ax_x_dist.twinx()
+    ax_xcum_dist.hist(x=column_corr, data=df, bins=40, cumulative=True, histtype='step',
+                      density=True, color='darkblue', align='mid')
+    ax_xcum_dist.set(xlim=(-1, 1))
+    ax_xcum_dist.tick_params(column_corr, colors='darkblue')
+    ax_xcum_dist.set_ylabel('cumulative', color='darkblue')
+    ax_xcum_dist.set(yticks=np.arange(0.0, 1.2, 0.2))
+
+    ax_y_dist.hist(x=column_pvalue, data=df, bins=200, orientation='horizontal',
+                   align='mid', color="lightblue")
+    ax_y_dist.set(xlabel='count')
+    ax_ycum_dist = ax_y_dist.twiny()
+    ax_ycum_dist.hist(x=column_pvalue, data=df, bins=200, cumulative=True, histtype='step',
+                      density=True, color='darkblue', align='mid', orientation='horizontal')
+    ax_ycum_dist.tick_params(column_pvalue, colors='darkblue')
+    ax_ycum_dist.set_xlabel('cumulative', color='darkblue')
+    ax_ycum_dist.set(xticks=np.arange(0.0, 1.2, 0.2), ylim=(0, 0.1))
+
+    plt.setp(ax_x_dist.get_xticklabels(), visible=False)
+    plt.setp(ax_y_dist.get_yticklabels(), visible=False)
+
+    plt.setp(ax_main.get_xticklabels(), rotation=90)
+    plt.setp(ax_y_dist.get_xticklabels(), rotation=90)
+
+    fig.subplots_adjust(top=0.85, right=0.85)
+    cbar_ax = fig.add_axes([0.87, 0.15, 0.03, 0.4])
+
+    cb = plt.colorbar(hb, cax=cbar_ax)
+    cb.set_label('counts')
+
+    return plt
+
+
+def plot_annotations(column_ppm_error, column_adducts, df):
 
     fig = plt.figure()
+    fig.set_size_inches(8.27, 11.69)
 
-    gs = gridspec.GridSpec(2, 2, height_ratios=[1, 5])
+    gs = gridspec.GridSpec(5, 2, height_ratios=[1, 1, 5, 1, 1])
 
-    ax_box = plt.subplot(gs[0])
-    ax_hist = plt.subplot(gs[2], sharex=ax_box)
-    ax_count = plt.subplot(gs[3])
-    #ax = plt.subplot(gs[1])
+    ax_box = plt.subplot(gs[2])
+    ax_hist = plt.subplot(gs[4], sharex=ax_box)
+    ax_count = plt.subplot(gs[5])
+    # ax = plt.subplot(gs[1])
 
     ppm_errors = df[column_ppm_error].dropna()
 
     sns.boxplot(ppm_errors, ax=ax_box)
     sns.distplot(ppm_errors, ax=ax_hist)
 
+    ax_hist.grid(False)
+    ax_box.grid(False)
+    ax_hist.grid(False)
+
     std = ppm_errors.std()
     mean = ppm_errors.mean()
     median = ppm_errors.median()
@@ -49,7 +129,7 @@ def report(df, column_ppm_error, column_adducts, fn_pdf):
 
     # Remove x axis name for the boxplot
     ax_box.set(xlabel="")
-    #ax_box.set_xticks([])
+    # ax_box.set_xticks([])
     ax_box.set_title("Q1={}; median={}; Q3={}".format(round(Q1, 2), round(median, 2), round(Q3, 2)))
 
     ax_hist.set_title("mean={}; std={}".format(round(mean, 2), round(std, 2)))
@@ -60,6 +140,33 @@ def report(df, column_ppm_error, column_adducts, fn_pdf):
 
     plt.setp(ax_box.get_xticklabels(), visible=False)
 
-    fig.suptitle('Summary - BEAMS', fontsize=20)
-    fig.set_size_inches(11.69, 8.27)
-    fig.savefig(fn_pdf, format="pdf")
+    return plt
+
+
+def report(db, pdf_out, column_corr, column_pvalue, column_ppm_error, column_adducts):
+
+    with PdfPages(pdf_out) as pdf:
+
+        conn = sqlite3.connect(db)
+        cursor = conn.cursor()
+        cursor.execute("""SELECT name FROM sqlite_master WHERE type='table';""")
+
+        for table in cursor.fetchall():
+            if str(table[0]) == "groups":
+
+                df = pd.read_sql_query("SELECT {}, {} FROM groups".format(column_corr, column_pvalue), conn)
+
+                plt = plot_correlations(column_corr, column_pvalue, df)
+                plt.suptitle('Summary - BEAMS\n\n\nGrouping features', fontsize=20)
+                pdf.savefig(dpi=300)
+                plt.close()
+
+            elif table[0][0:10] == "compounds_":
+
+                df = pd.read_sql_query("SELECT {}, {} FROM {}".format(column_ppm_error, column_adducts, table[0]), conn)
+
+                plt = plot_annotations("ppm_error", "adduct", df)
+                plt.suptitle('Annotation\n\n{}'.format(table[0].replace("compounds_", "")), fontsize=20)
+                pdf.savefig(dpi=300)
+                plt.close()
+        conn.close()

conda/meta.yaml

@@ -5,30 +5,25 @@
 import os
 import numpy as np
 import pandas as pd
-from tests.utils import to_test_results
+from tests.utils import to_test_data, to_test_results
 from beams.plots import report
 
 
 class PlotsTestCase(unittest.TestCase):
 
     def test_report(self):
 
-        np.random.seed(0)
-        n = 1000
-        mu, sigma = 0, 0.1 # mean and standard deviation
-        s = np.random.normal(mu, sigma, n)
-
-        lib = ["[M+H]+", "[M+Na]+", "[M+K]+"]
-        adducts = [lib[i] for i in np.random.randint(3, size=n)]
-
-        df = pd.DataFrame({'ppm_error': s, "adduct": adducts})
-
-        report(df, "ppm_error", "adduct", to_test_results("test_report_01.pdf"))
-
+        report(to_test_data("results_annotation.sqlite"), to_test_results("test_report_01.pdf"),
+               "r_value", "p_value", "ppm_error", "adduct")
         statinfo = os.stat(to_test_results("test_report_01.pdf"))
-        # print statinfo.st_size
-        self.assertTrue(statinfo.st_size > 16000)
-
+        #print(statinfo.st_size)
+        self.assertTrue(statinfo.st_size > 700)
+
+        report(to_test_data("results_pearson_all.sqlite"),  to_test_results("test_report_02.pdf"),
+               "r_value", "p_value", "ppm_error", "adduct")
+        statinfo = os.stat(to_test_results("test_report_02.pdf"))
+        print(statinfo.st_size)
+        #self.assertTrue(statinfo.st_size > 25000)
 
 if __name__ == '__main__':
     unittest.main()