Seurat VlnPlot presenting expression of multiple genes in a single cluster
To do so one workaround it to have your data in "long format" and then use the column that holds the "gene names" as the x variable while plotting.

You can use FetchData() to extract data from a Seurat object. VlnPlot's default is the data slot (of the active assay if using Seurat v3 I suppose). And you can specify which cells and genes to retrieve.

selected_cells &lt;- names(panc8$celltype[panc8$celltype == "gamma"])
data &lt;- FetchData(panc8,
                    vars = c('FZD9','CTNNB1','APC'),
                    cells = selected_cells ,
                    slot = "data")
&gt; head(data)
           FZD9   CTNNB1      APC
D101_5        0 0.000000 0.000000
D101_43       0 2.007853 1.001958
D101_93       0 0.000000 0.000000


Melt() would transform your data into "long format". By not specifying any arguments, we push for all of the info in the three variables to be gathered in two columns:

long_data &lt;- melt(data)
No id variables; using all as measure variables
&gt; head(long_data)
  variable value
1     FZD9     0
2     FZD9     0
3     FZD9     0

&gt; tail(long_data)
     variable value
1870      APC     0
1871      APC     0
1872      APC     0


ggplot2 is used to add "violin" and "jitter" layers. You can customize the output to look (exactly) like the VlnPlot() output.

ggplot(long_data,
       aes(x = variable, y = value)) +
  geom_violin() +
  geom_jitter(size = 0.1)




And here is the same graph generated by VlnPlot():



There are no "violins" as the counts are almost entirely zeros. And see how different the ranges are in the y axes of the VlnPlot. So the "tweak" I have presented here would only work for genes that are expressed at similar levels / similar ranges.