docs: def volcano(df, a, b=None,
              fdr=0.05,
              figsize=(8,10),
              show_numbers=True,
              threshold=2,
              minimum_sample_n=0,
              estimate_qvalues=False,
              labels_from=None,
              labels_for=None,
              title=None,
              label_format=None,
              markersize=64,
              s0=0.00001,
              draw_fdr=True,
              is_log2=False,
              fillna=None,
              label_sig_only=True,
              ax=None,
              xlim=None,
              ylim=None,
              fc='grey',
              fc_sig='blue',
              fc_sigr='red',
              optimize_label_iter=OPTIMIZE_LABEL_ITER_DEFAULT
              ):
      """
      Volcano plot of two sample groups showing t-test p value vs. log2(fc).
  
      Generates a volcano plot for two sample groups, selected from `df` using `a` and `b` indexers. The mean of
      each group is calculated along the y-axis (per protein) and used to generate a log2 ratio. If a log2-transformed
      dataset is supplied set `islog2=True` (a warning will be given when negative values are present).
  
      If `a` and `b` are provided a two-sample independent t-test is performed between each group. If `minimum_sample_n` is supplied, any values (proteins)
      without this number of samples will be dropped from the analysis.
  
      If only `a` is provided, a one-sample t-test is performend.
  
      Individual data points can be labelled in the resulting plot by passing `labels_from` with a index name, and `labels_for`
      with a list of matching values for which to plot labels.
  
  
      :param df: Pandas `dataframe`
      :param a: `tuple` or `str` indexer for group A
      :param b: `tuple` or `str` indexer for group B
      :param fdr: `float` false discovery rate cut-off
      :param threshold: `float` log2(fc) ratio cut -off
      :param minimum_sample_n: `int` minimum sample for t-test
      :param estimate_qvalues: `bool` estimate Q values (adjusted P)
      :param labels_from: `str` or `int` index level to get labels from
      :param labels_for: `list` of `str` matching labels to show
      :param title: `str` title for plot
      :param markersize: `int` size of markers
      :param s0: `float` smoothing factor between fdr/fc cutoff
      :param draw_fdr: `bool` draw the fdr/fc curve
      :param is_log2: `bool` is the data log2 transformed already?
      :param fillna: `float` fill NaN values with value (default: 0)
      :param label_sig_only: `bool` only label significant values
      :param ax: matplotlib `axis` on which to draw
      :param fc: `str` hex or matplotlib color code, default color of points
      :param optimize_label_iter: `50` (default) number of iterations to attempt to optimize label positions.
                                  Increase for clearer positions (longer processing time).
      :return:
      """
      df = df.copy()
      
      if np.any(df.values[~pd.isnull(df.values)] < 0) and not is_log2:
          warnings.warn("Input data has negative values. If data is log2 transformed, set is_log2=True.")
  
      if fillna is not None:
          df = df.fillna(fillna)
          
      if labels_from is None:
          labels_from = list(df.index.names)
  
      if b is not None:
          # Calculate ratio between two groups
          g1, g2 = df[a].values, df[b].values
  
          if is_log2:
              dr = np.nanmean(g2, axis=1) - np.nanmean(g1, axis=1)
          else:
              dr = np.log2(np.nanmean(g2, axis=1) / np.nanmean(g1, axis=1))
  
          ginv = ( (~np.isnan(g1)).sum(axis=1) < minimum_sample_n ) | ((~np.isnan(g2)).sum(axis=1) < minimum_sample_n)
  
          g1 = np.ma.masked_where(np.isnan(g1), g1)
          g2 = np.ma.masked_where(np.isnan(g2), g2)
  
          # Calculate the p value between two groups (t-test)
          t, p = sp.stats.mstats.ttest_ind(g1.T, g2.T)
  
      else:
          g1 = df[a].values
          dr = np.nanmean(g1, axis=1)
  
          ginv = (~np.isnan(g1)).sum(axis=1) < minimum_sample_n
  
          # Calculate the p value one sample t
          g1 = np.ma.masked_where(np.isnan(g1), g1)
          t, p = sp.stats.mstats.ttest_1samp(g1.T, popmean=0)
  
      p = np.array(p) # Unmask
  
      # Set p values to nan where not >= minimum_sample_n values
      p[ ginv ] = np.nan
  
      if ax is None:
          fig = plt.figure(figsize=figsize)
  
          ax = fig.add_subplot(1,1,1)
  
      if estimate_qvalues:
          p[~np.isnan(p)] = qvalues(p[~np.isnan(p)])
          ax.set_ylabel('-log$_{10}$(Q)')
      else:
          ax.set_ylabel('-log$_{10}$(p)')
  
  
      # There are values below the fdr
      s0x, s0y, s0fn = calculate_s0_curve(s0, fdr, min(fdr/2, np.nanmin(p)), np.log2(threshold), np.nanmax(np.abs(dr)), curve_interval=0.001)
  
      if draw_fdr is True:
          ax.plot(s0x, -np.log10(s0y), fc_sigr, lw=1 )
          ax.plot(-s0x, -np.log10(s0y), fc_sigr, lw=1 )
  
      # Select data based on s0 curve
      _FILTER_IN = []
  
      for x, y in zip(dr, p):
          x = np.abs(x)
          spy = s0fn(x)
  
          if len(s0x) == 0 or x < np.min(s0x):
              _FILTER_IN.append(False)
              continue
  
          if y <= spy:
              _FILTER_IN.append(True)
          else:
              _FILTER_IN.append(False)
  
      _FILTER_IN = np.array(_FILTER_IN)
      _FILTER_OUT = ~ _FILTER_IN
      
  
      def scatter(ax, f, c, alpha=0.5):
          # lw = [float(l[0][-1])/5 for l in df[ f].index.values]
          
          if type(markersize) == str:
              # Use as a key vs. index value in this levels
              s = df.index.get_level_values(markersize)
          elif callable(markersize):
              s = np.array([markersize(c) for c in df.index.values])
          else:
              s = np.ones((df.shape[0],))*markersize
          
          
          ax.scatter(dr[f], -np.log10(p[f]), c=c, s=s[f], linewidths=0, alpha=0.5)
      
      _FILTER_OUT1 = _FILTER_OUT & ~np.isnan(p) & (np.array(p) > fdr)
      scatter(ax, _FILTER_OUT1, fc, alpha=0.3)
  
      _FILTER_OUT2 = _FILTER_OUT & (np.array(p) <= fdr)
      scatter(ax, _FILTER_OUT2, fc_sig, alpha=0.3)
  
      if labels_for:
          texts = []
          idxs = get_index_list( df.index.names, labels_from )
          for shown, label, x, y in zip( _FILTER_IN , df.index.values, dr, -np.log10(p)):
              
              if shown or not label_sig_only:
                  label = build_combined_label( label, idxs, label_format=label_format)
                  
                  if labels_for == True or any([l in label for l in labels_for]):
                      t = ax.text(x, y, label , ha='center', va='center', rotation_mode='anchor', bbox=dict(boxstyle='round,pad=0.3', fc='#ffffff', ec='none', alpha=0.4))
                      texts.append(t)
  
          # Adjust spacing/positioning of labels if required.
          if texts and optimize_label_iter:
              adjust_text(
                  texts,
                  lim=optimize_label_iter,
                  arrowprops=dict(arrowstyle='->', color='red')
              )
  
      scatter(ax, _FILTER_IN, fc_sigr)
      
  
      ax.set_xlabel('log$_2$ ratio')
  
      # Centre the plot horizontally
      if xlim is None:
          xmin, xmax = ax.get_xlim()
          xlim = np.max(np.abs([xmin, xmax]))
      ax.set_xlim((-xlim, xlim))
  
      if ylim is None:
          _, ylim = ax.get_ylim()
      ax.set_ylim((0, ylim))
  
      if title:
          ax.set_title(title)
  
      if show_numbers:
          # Annotate axes with the numbers of points in each category (up & down):
          # filtered (red), Sig-filtered (blue), nonsig (grey)
          dr_up, dr_down = dr > 0, dr < 0
          grey_up, blue_up, red_up = np.sum(_FILTER_OUT1 & dr_up), np.sum(_FILTER_OUT2 & dr_up), np.sum(_FILTER_IN & dr_up)
          grey_do, blue_do, red_do = np.sum(_FILTER_OUT1 & dr_down), np.sum(_FILTER_OUT2 & dr_down), np.sum(_FILTER_IN & dr_down)
  
          ax.text(0.95, 0.95, '%d' % red_up, horizontalalignment='right', transform=ax.transAxes, color=fc_sigr, fontsize=14)
          ax.text(0.95, 0.90, '%d' % blue_up, horizontalalignment='right', transform=ax.transAxes, color=fc_sig, fontsize=14)
          ax.text(0.95, 0.05, '%d' % grey_up, horizontalalignment='right', transform=ax.transAxes, color=fc, fontsize=14)
  
          ax.text(0.05, 0.95, '%d' % red_do, horizontalalignment='left', transform=ax.transAxes, color=fc_sigr, fontsize=14)
          ax.text(0.05, 0.90, '%d' % blue_do, horizontalalignment='left', transform=ax.transAxes, color=fc_sig, fontsize=14)
          ax.text(0.05, 0.05, '%d' % grey_do, horizontalalignment='left', transform=ax.transAxes, color=fc, fontsize=14)
  
      return ax, p, dr, _FILTER_IN