sc.pl.violin doesn't work with `groupby` and `ax` passed together

• I have checked that this issue has not already been reported.
• I have confirmed this bug exists on the latest version of scanpy.
• (optional) I have confirmed this bug exists on the master branch of scanpy.

Minimal code sample (that we can copy&paste without having any data)

Option 1: group with two keys, passing two axes

import scanpy as sc
import matplotlib.pyplot as plt
adata = sc.datasets.pbmc3k()
adata.obs['group'] = adata.obs.index.to_series().str.startswith("A").astype(str)
fig, axes = plt.subplots(1, 2)
sc.pl.violin(adata2, keys=['CD8A', 'CD8B'], groupby="group", ax=axes)

---------------------------------------------------------------------------
AttributeError                            Traceback (most recent call last)
Input In [51], in <module>
      4 adata2.obs['group'] = adata2.obs.index.to_series().str.startswith("A").astype(str)
      5 fig, axes = plt.subplots(1, 2)
----> 6 sc.pl.violin(adata2, keys=['CD8A', 'CD8B'], groupby="group", ax=axes)

File /opt/conda/envs/analysis/lib/python3.8/site-packages/scanpy/plotting/_anndata.py:835, in violin(adata, keys, groupby, log, use_raw, stripplot, jitter, size, layer, scale, order, multi_panel, xlabel, ylabel, rotation, show, save, ax, **kwds)
    833     axs = [ax]
    834 for ax, y, ylab in zip(axs, ys, ylabel):
--> 835     ax = sns.violinplot(
    836         x=x,
    837         y=y,
    838         data=obs_tidy,
    839         order=order,
    840         orient='vertical',
    841         scale=scale,
    842         ax=ax,
    843         **kwds,
    844     )
    845     if stripplot:
    846         ax = sns.stripplot(
    847             x=x,
    848             y=y,
   (...)
    854             ax=ax,
    855         )

File /opt/conda/envs/analysis/lib/python3.8/site-packages/seaborn/_decorators.py:46, in _deprecate_positional_args.<locals>.inner_f(*args, **kwargs)
     36     warnings.warn(
     37         "Pass the following variable{} as {}keyword arg{}: {}. "
     38         "From version 0.12, the only valid positional argument "
   (...)
     43         FutureWarning
     44     )
     45 kwargs.update({k: arg for k, arg in zip(sig.parameters, args)})
---> 46 return f(**kwargs)

File /opt/conda/envs/analysis/lib/python3.8/site-packages/seaborn/categorical.py:2408, in violinplot(x, y, hue, data, order, hue_order, bw, cut, scale, scale_hue, gridsize, width, inner, split, dodge, orient, linewidth, color, palette, saturation, ax, **kwargs)
   2405 if ax is None:
   2406     ax = plt.gca()
-> 2408 plotter.plot(ax)
   2409 return ax

File /opt/conda/envs/analysis/lib/python3.8/site-packages/seaborn/categorical.py:1043, in _ViolinPlotter.plot(self, ax)
   1041 def plot(self, ax):
   1042     """Make the violin plot."""
-> 1043     self.draw_violins(ax)
   1044     self.annotate_axes(ax)
   1045     if self.orient == "h":

File /opt/conda/envs/analysis/lib/python3.8/site-packages/seaborn/categorical.py:761, in _ViolinPlotter.draw_violins(self, ax)
    759 def draw_violins(self, ax):
    760     """Draw the violins onto `ax`."""
--> 761     fill_func = ax.fill_betweenx if self.orient == "v" else ax.fill_between
    762     for i, group_data in enumerate(self.plot_data):
    764         kws = dict(edgecolor=self.gray, linewidth=self.linewidth)

AttributeError: 'numpy.ndarray' object has no attribute 'fill_betweenx'

Option 2: group with two keys, passing first of two axes

import scanpy as sc
import matplotlib.pyplot as plt
adata = sc.datasets.pbmc3k()
adata.obs['group'] = adata.obs.index.to_series().str.startswith("A").astype(str)
fig, axes = plt.subplots(1, 2)
sc.pl.violin(adata2, keys=['CD8A', 'CD8B'], groupby="group", ax=axes[0])

No traceback, but the second axis is simply not plotted.

Option 3: group with two keys, passing one axis

import scanpy as sc
import matplotlib.pyplot as plt
adata = sc.datasets.pbmc3k()
adata.obs['group'] = adata.obs.index.to_series().str.startswith("A").astype(str)
fig, ax = plt.subplots()
sc.pl.violin(adata, keys=['CD8A', 'CD8B'], groupby="group", ax=ax)

No traceback, even though this should error. Plots just the first of the two keys.

Versions

WARNING: If you miss a compact list, please try `print_header`!
-----
anndata     0.7.8
scanpy      1.8.2
sinfo       0.3.1
-----
PIL                         8.4.0
anndata                     0.7.8
asttokens                   NA
attr                        21.2.0
backcall                    0.2.0
cffi                        1.15.0
colorama                    0.4.4
cycler                      0.10.0
cython_runtime              NA
dateutil                    2.8.0
debugpy                     1.5.1
decorator                   5.1.1
defusedxml                  0.7.1
django                      4.0
executing                   0.8.2
google                      NA
h5py                        2.10.0
idna                        3.1
igraph                      0.9.8
importlib_resources         NA
ipykernel                   6.7.0
ipython_genutils            0.2.0
ipywidgets                  7.6.5
jedi                        0.18.1
jinja2                      3.0.3
joblib                      1.1.0
jsonschema                  4.4.0
jupyter_server              1.13.3
kiwisolver                  1.3.2
leidenalg                   0.8.8
llvmlite                    0.37.0
markupsafe                  2.0.1
matplotlib                  3.5.1
matplotlib_inline           NA
mpl_toolkits                NA
natsort                     8.0.2
nbformat                    5.1.3
numba                       0.54.1
numexpr                     2.8.0
numpy                       1.19.5
packaging                   21.3
pandas                      1.1.5
parso                       0.8.3
pexpect                     4.8.0
pickleshare                 0.7.5
pkg_resources               NA
prometheus_client           NA
prompt_toolkit              3.0.24
psutil                      5.8.0
ptyprocess                  0.7.0
pure_eval                   0.2.1
pvectorc                    NA
pydev_ipython               NA
pydevconsole                NA
pydevd                      2.6.0
pydevd_concurrency_analyser NA
pydevd_file_utils           NA
pydevd_plugins              NA
pydevd_tracing              NA
pygments                    2.10.0
pyparsing                   3.0.6
pyrsistent                  NA
pytz                        2021.3
scanpy                      1.8.2
scipy                       1.5.3
scprep                      1.1.0
seaborn                     0.11.2
send2trash                  NA
setuptools                  58.0.4
setuptools_scm              NA
sinfo                       0.3.1
six                         1.16.0
sklearn                     0.24.2
sphinxcontrib               NA
stack_data                  0.1.4
statsmodels                 0.13.1
tables                      3.6.1
terminado                   0.12.1
texttable                   1.6.4
tornado                     6.1
tqdm                        4.62.3
traitlets                   5.1.1
typing_extensions           NA
wcwidth                     0.2.5
yaml                        6.0
zipp                        NA
zmq                         22.3.0
-----
IPython             8.0.0
jupyter_client      6.1.12
jupyter_core        4.9.1
jupyterlab          3.2.8
notebook            6.4.7
-----
Python 3.8.12 | packaged by conda-forge | (default, Oct 12 2021, 21:59:51) [GCC 9.4.0]
Linux-5.4.0-1064-gcp-x86_64-with-glibc2.10
16 logical CPU cores
-----
Session information updated at 2022-02-10 16:38

Start with scanpy/plotting/_anndata.py, the file shown in the traceback, and run the minimal reproduction for sc.pl.violin with groupby and ax. Check how the three ax cases are handled; done means multiple keys plot on the intended axes, invalid combinations fail clearly, and regression coverage is added.