src/fig_utils/: modules to add statistical annoation on figure

src/fig_utils/__init__.py

+#!/usr/bin/env python3
+
+# -*- coding: UTF-8 -*-
+
+"""
+Description:
+"""
\ No newline at end of file

src/fig_utils/stat_annot.py

+#!/usr/bin/env python3
+
+# -*- coding: UTF-8 -*-
+
+"""
+Description:
+    Add annotation on figure.
+    Script adapted from
+    github.com/webermarcolivier/statannot/blob/master/statannot/statannot.py.
+"""
+
+
+import matplotlib.pyplot as plt
+from matplotlib import lines
+import matplotlib.transforms as mtransforms
+from matplotlib.font_manager import FontProperties
+import numpy as np
+import seaborn as sns
+from seaborn.utils import remove_na
+
+
+def add_stat_annotation(ax, data=None, x=None, y=None, hue=None, order=None,
+                        hue_order=None, box_pair_list=None, loc='inside',
+                        use_fixed_offset=False,
+                        liney_offset2box_axes_coord=None,
+                        liney_offset_axes_coord=None,
+                        line_height_axes_coord=0.02,
+                        text_y_offset_points=1, color='0.2', linewidth=1.5,
+                        fontsize='medium', verbose=1):
+    """
+    User should use the same argument for the data, x, y, hue, order, \
+    hue_order as the seaborn boxplot function.
+    box_pair_list can be of either form:
+    boxplot: [(cat1, cat2, pval), (cat3, cat4, pval)]
+    """
+
+    def find_x_position_box(box_plotter, box_name):
+        """
+        box_name can be either a name "cat" or a tuple ("cat", "hue")
+        """
+        if box_plotter.plot_hues is None:
+            cat = box_name
+            hue_offset = 0
+        else:
+            cat = box_name[0]
+            chue = box_name[1]
+            hue_offset = box_plotter.hue_offsets[
+                box_plotter.hue_names.index(chue)]
+
+        group_pos = box_plotter.group_names.index(cat)
+        box_pos = group_pos + hue_offset
+        return box_pos
+
+    def get_box_data(cbox_plotter, box_name):
+        """
+        box_name can be either a name "cat" or a tuple ("cat", "hue")
+        Here we really have to duplicate seaborn code, because there is not
+        direct access to the
+        box_data in the BoxPlotter class.
+        """
+        cat = box_name
+        if cbox_plotter.plot_hues is None:
+            box_max_l = []
+            for cat in cbox_plotter.group_names:
+                i = cbox_plotter.group_names.index(cat)
+                group_data = cbox_plotter.plot_data[i]
+                box_data = remove_na(group_data)
+                box_max_l.append(np.mean(box_data))
+            return max(box_max_l)
+        else:
+            i = cbox_plotter.group_names.index(cat[0])
+            group_data = cbox_plotter.plot_data[i]
+            box_data = []
+            for chue in cbox_plotter.hue_names:
+                hue_mask = cbox_plotter.plot_hues[i] == chue
+                box_data.append(np.max(remove_na(group_data[hue_mask])))
+            return max(box_data)
+
+    fig = plt.gcf()
+
+    valid_list = ['inside', 'outside']
+    if loc not in valid_list:
+        raise ValueError(f"loc value should be one of the following: "
+                         f"{', '.join(valid_list)}.")
+
+    # Create the same BoxPlotter object as seaborn's boxplot
+    box_plotter = sns.categorical._BoxPlotter(x, y, hue, data, order,
+                                              hue_order, orient=None, width=.8,
+                                              color=None, palette=None,
+                                              saturation=.75, dodge=True,
+                                              fliersize=5, linewidth=None)
+    ylim = ax.get_ylim()
+    y_range = ylim[1] - ylim[0]
+
+    if liney_offset_axes_coord is None:
+        if loc == 'inside':
+            liney_offset_axes_coord = 0.005
+            if liney_offset2box_axes_coord is None:
+                liney_offset2box_axes_coord = 0.1
+        elif loc == 'outside':
+            liney_offset_axes_coord = 0.03
+            liney_offset2box_axes_coord = liney_offset_axes_coord
+    else:
+        if loc == 'inside':
+            if liney_offset2box_axes_coord is None:
+                liney_offset2box_axes_coord = 0.06
+        elif loc == 'outside':
+            liney_offset2box_axes_coord = liney_offset_axes_coord
+    y_offset = liney_offset_axes_coord * y_range
+    y_offset_to_box = liney_offset2box_axes_coord * y_range
+
+    y_stack = []
+    ann_list = []
+    for box1, box2, pval in box_pair_list:
+
+        group_names = box_plotter.group_names
+        cat1 = box1
+        cat2 = box2
+        if isinstance(cat1, tuple):
+            hue_names = box_plotter.hue_names
+            valid = cat1[0] in group_names and cat2[0] in group_names and \
+                cat1[1] in hue_names and cat2[1] in hue_names
+        else:
+            valid = cat1 in group_names and cat2 in group_names
+
+        if valid:
+            # Get position of boxes
+            x1 = find_x_position_box(box_plotter, box1)
+            x2 = find_x_position_box(box_plotter, box2)
+            box_data1 = get_box_data(box_plotter, box1)
+            box_data2 = get_box_data(box_plotter, box2)
+            ymax1 = box_data1
+            ymax2 = box_data2
+
+            if pval > 1e-16:
+                text = "p = {:.2e}".format(pval)
+            else:
+                text = "p < 1e-16"
+
+            if loc == 'inside':
+                y_ref = max(ymax1, ymax2)
+            else:
+                y_ref = ylim[1]
+            if len(y_stack) > 0:
+                y_ref2 = max(y_ref, max(y_stack))
+            else:
+                y_ref2 = y_ref
+
+            if len(y_stack) == 0:
+                y = y_ref2 + y_offset_to_box
+            else:
+                y = y_ref2 + y_offset
+            h = line_height_axes_coord * y_range
+            line_x, line_y = [x1, x1, x2, x2], [y, y + h, y + h, y]
+            if loc == 'inside':
+                ax.plot(line_x, line_y, lw=linewidth, c=color)
+            elif loc == 'outside':
+                line = lines.Line2D(line_x, line_y, lw=linewidth, c=color,
+                                    transform=ax.transData)
+                line.set_clip_on(False)
+                ax.add_line(line)
+
+            if text is not None:
+                ann = ax.annotate(text, xy=(np.mean([x1, x2]), y + h),
+                                  xytext=(0, text_y_offset_points),
+                                  textcoords='offset points',
+                                  xycoords='data', ha='center', va='bottom',
+                                  fontsize=fontsize, clip_on=False,
+                                  annotation_clip=False)
+                ann_list.append(ann)
+
+            new_max_ylim = 1.1*(y + h)
+            if new_max_ylim > ylim[1]:
+                ax.set_ylim((ylim[0], 1.1*(y + h)))
+
+            if text is not None:
+                plt.draw()
+                y_top_annot = None
+                got_matplotlib_error = False
+                if not use_fixed_offset:
+                    try:
+                        bbox = ann.get_window_extent()
+                        bbox_data = bbox.transformed(ax.transData.inverted())
+                        y_top_annot = bbox_data.ymax
+                    except RuntimeError:
+                        got_matplotlib_error = True
+
+                if use_fixed_offset or got_matplotlib_error:
+                    if verbose >= 1:
+                        print("Warning: cannot get the text bounding box. "
+                              "Falling back to a fixed y offset. Layout may "
+                              "be not optimal.")
+                    #  We will apply a fixed offset in points,
+                    #  based on the font size of the annotation.
+                    fontsize_points = FontProperties(size='medium'
+                                                     ).get_size_in_points()
+                    offset_trans = mtransforms.offset_copy(
+                        ax.transData, fig=fig, x=0,
+                        y=1.0 * fontsize_points + text_y_offset_points,
+                        units='points')
+                    y_top_display = offset_trans.transform((0, y + h))
+                    y_top_annot = ax.transData.inverted().\
+                        transform(y_top_display)[1]
+            else:
+                y_top_annot = y + h
+
+            y_stack.append(y_top_annot)
+        else:
+            raise ValueError("box_pair_list contains an unvalid box pair.")
+
+    y_stack_max = max(y_stack)
+    if loc == 'inside':
+        if ylim[1] < 1.03 * y_stack_max:
+            ax.set_ylim((ylim[0], 1.03 * y_stack_max))
+    elif loc == 'outside':
+        ax.set_ylim((ylim[0], ylim[1]))
+    return ax