How to Make Violin Plots of RNA-seq Expression Data

By the end of this walkthrough you will have a violin plot you can defend. Expression values sit on a stated scale of log2 CPM or log-normalized counts. The kernel density uses a bandwidth you chose on purpose, zeros are handled explicitly rather than smoothed away, and per-group sample counts are printed on the axis. The tooling is either Python 3.10+ with scanpy, anndata, seaborn, matplotlib, pandas, and numpy, or R with DESeq2, ggplot2, and, for single-cell work, Seurat.

The inputs differ by assay. For bulk RNA-seq you need a raw count matrix (genes × samples, integers) together with a sample metadata table. For single cell you need a filtered feature-barcode matrix (matrix.mtx.gz, features.tsv.gz, barcodes.tsv.gz) or an .h5ad file. If you have your own RNA-seq from a profiling service, you are almost certainly in the bulk case, holding either a salmon or kallisto quantification directory or a gene-level counts TSV. Start there.

It helps to be clear about what the figure is. A violin plot is a mirrored kernel density estimate of a distribution, usually annotated with quartile marks or an inner box. It answers one question well: what is the shape of expression for this gene across these cells or samples, and is that shape unimodal. It answers “is this gene differentially expressed” badly. It shows no model, no dispersion estimate, and no multiple-testing correction. Treat it as a sanity check on a statistical result rather than as the result itself.

1. Get expression onto a scale that makes density estimation meaningful

Everything downstream depends on the scale you plot, so this is the first decision to make deliberately. Raw counts are the wrong choice. They are dominated by library size and their variance scales with the mean, which means a kernel density estimate over raw counts is largely a picture of sequencing depth. Decide what you are plotting before you plot it.

For bulk RNA-seq, use DESeq2’s variance-stabilizing transform or log2 CPM. The variance-stabilizing transform (VST) is the better default when you have at least 10 samples, because it flattens the mean-variance relationship, and that flattening is what makes comparing genes on a shared axis defensible. In R:

library(DESeq2)
cts <- as.matrix(read.delim("counts.tsv", row.names = 1))
coldata <- read.delim("samples.tsv", row.names = 1)
stopifnot(all(colnames(cts) == rownames(coldata)))

dds <- DESeqDataSetFromMatrix(cts, coldata, design = ~ condition)
dds <- dds[rowSums(counts(dds) >= 10) >= 3, ]   # drop near-zero genes
vsd <- vst(dds, blind = FALSE)                   # blind=TRUE for pure QC
expr <- assay(vsd)                               # log2-ish scale

Smaller cohorts need a different treatment. With fewer than about 10 samples, vst() falls back on a parametric fit that can be unstable, so use rlog(dds, blind = FALSE) instead, or plain log2(cpm + 1) via edgeR’s cpm(dds_counts, log = TRUE, prior.count = 2). The prior.count = 2 matters more than it looks: with prior.count = 0.5 or less, low-count genes produce a long left tail of large negative values, and the violin grows a spurious lower lobe that is pure shrinkage artifact.

For single-cell data, plot the log-normalized layer rather than raw UMIs or scaled (z-scored) values. Scaled values are centered per gene, so the violin sits at zero and the y-axis stops meaning anything biological.

import scanpy as sc
adata = sc.read_10x_mtx("filtered_feature_bc_matrix/", var_names="gene_symbols", cache=True)
adata.layers["counts"] = adata.X.copy()
sc.pp.normalize_total(adata, target_sum=1e4)
sc.pp.log1p(adata)
adata.raw = adata           # keeps log-normalized values for plotting

Plot from adata.raw or from an explicit layer. If you later run sc.pp.scale, sc.pl.violin will silently plot the scaled values unless you pass use_raw=True, and that is the single most common way people publish a meaningless y-axis.

2. Make the QC violins first

Before you plot any individual gene, look at the three quality-control metrics as violins. This is the plot the Seurat tutorial shows. The black dots people ask about on Stack Exchange are the individual cells, jittered horizontally, with the violin drawn as their density.

adata.var["mt"] = adata.var_names.str.startswith(("MT-", "mt-"))
sc.pp.calculate_qc_metrics(adata, qc_vars=["mt"], percent_top=None,
                           log1p=False, inplace=True)

sc.pl.violin(adata, ["n_genes_by_counts", "total_counts", "pct_counts_mt"],
             jitter=0.4, multi_panel=True, stripplot=True, size=1)

Each panel carries a specific signal. A sharp left edge on total_counts means your cell-calling threshold cut there rather than reflecting biology. A second lobe at roughly twice the main mode indicates doublets when it appears in both total_counts and n_genes_by_counts. Run scrublet or DoubletFinder before going further. A heavy right tail on pct_counts_mt past roughly 15-20% in most tissues points to stressed or dying cells. Choose your thresholds from the plot, write them down, and apply them once:

adata = adata[(adata.obs.n_genes_by_counts > 500)
              & (adata.obs.pct_counts_mt < 15)].copy()

Resist the temptation to pick thresholds from a published table of standard values. Nuclei preparations legitimately sit at 1-2% mitochondrial reads, and cardiomyocytes legitimately sit high. The violin tells you what your own sample looks like.

3. Plot one gene across groups (bulk)

For a personal bulk dataset, this is usually the most useful figure. It shows a single gene and its expression across timepoints or conditions, with every sample visible as a point.

import pandas as pd, seaborn as sns, matplotlib.pyplot as plt

expr = pd.read_csv("vst.tsv", sep="\t", index_col=0)      # genes x samples
meta = pd.read_csv("samples.tsv", sep="\t", index_col=0)

gene = "HMOX1"
df = meta.join(expr.loc[gene].rename("expr"))

fig, ax = plt.subplots(figsize=(4.2, 3.4))
sns.violinplot(data=df, x="condition", y="expr", hue="condition",
               inner="quartile", cut=0, bw_adjust=1.0, density_norm="width",
               linewidth=1, legend=False, ax=ax)
sns.stripplot(data=df, x="condition", y="expr", color="black",
              size=3.5, jitter=0.18, alpha=0.9, ax=ax)

n = df.groupby("condition").size()
ax.set_xticks(range(len(n)))
ax.set_xticklabels([f"{k}\nn={v}" for k, v in n.items()])
ax.set_ylabel(f"{gene} (VST)")
fig.tight_layout()
fig.savefig("hmox1_violin.pdf")

Four arguments in that call carry all the weight, and each is worth understanding.

cut=0 truncates the density at the observed data range. The default, cut=2, extends the kernel density estimate two bandwidths past your minimum and maximum. That draws expression values you never measured, including negative ones on a log scale. Use cut=0 for expression data every time.

density_norm="width" (called scale="width" in older seaborn) gives every violin the same maximum width so that shapes are comparable across groups. Use "count" instead if group sizes differ substantially and you want the plot to show that difference. Avoid the default "area" if you plan to describe widths across groups, since a wide-but-flat group and a narrow-but-peaked group can have identical area.

bw_adjust scales the bandwidth, with seaborn using Scott’s rule by default. On small bulk cohorts of fewer than 8 samples per group, Scott’s rule oversmooths everything into a single fat blob. Try bw_adjust=0.6 and see whether structure appears; if structure shows up only at 0.4 and vanishes at 0.8, it is noise. Report the value you used in the figure caption.

inner="quartile" draws the 25th, 50th, and 75th percentile lines. With fewer than 10 samples those lines are noisy enough that better options exist. Try inner="stick", which draws one tick per observation, or inner=None combined with the stripplot. With n below about 6 per group, do not draw a violin at all. A kernel density estimate over five points is a picture of your bandwidth rather than your data, so use a dotplot or a beeswarm.

4. Plot one gene across clusters (single cell)

The single-cell analogue puts clusters on the x-axis and one gene on the y-axis.

sc.pl.violin(adata, keys=["SERPINE2"], groupby="leiden",
             use_raw=True, log=False, stripplot=False,
             rotation=90, cut=0, scale="width")

Turn stripplot off above a few thousand cells, because the dots merge into a solid black bar and hide the density you came to see. When you want many genes at once, stacked violins are far more readable than a grid of separate panels:

markers = ["PTPRC", "CD3D", "CD8A", "NKG7", "LYZ", "CD79A", "EPCAM", "PECAM1", "COL1A1"]
sc.pl.stacked_violin(adata, markers, groupby="leiden",
                     use_raw=True, swap_axes=False, standard_scale="var",
                     cmap="Blues", dendrogram=True)

standard_scale="var" rescales each gene to the range [0,1] across groups, which is what you want when the question is marker identity and the wrong choice if the reader will read absolute levels off the figure. If absolute levels matter, drop the argument and accept that high-expression genes will dominate the visual range.

This pattern puts clusters on the x-axis and one candidate gene on the y-axis, with marker violins supplying the identity evidence. It is how single-cell papers establish that a gene of interest is restricted to a particular cell compartment. The MYBL2 work in ovarian cancer uses exactly this approach to place expression in malignant epithelial cells rather than stroma or immune cells¹, and the SERPINE2 renal cell carcinoma study does the same before moving on to functional validation². The plot serves as the hand-off between clustering and the claim about which cell type carries the signal.

5. Handle the zeros

Single-cell violins have a structural problem worth addressing directly. Most cells have zero counts for most genes, so the kernel density estimate puts a large mass at exactly 0 and then smooths that mass into a lobe spreading below zero, since the estimator does not know your data is bounded. Setting cut=0 fixes the below-zero part. The mass at zero is real information, and you should quantify it rather than let the violin’s shape stand in for it.

import numpy as np
X = adata.raw[:, "SERPINE2"].X
v = np.asarray(X.todense()).ravel() if hasattr(X, "todense") else np.ravel(X)
for g, idx in adata.obs.groupby("leiden").indices.items():
    x = v[idx]
    print(f"{g}\tn={len(x)}\tdetected={np.mean(x > 0):.2%}\t"
          f"mean_expressing={x[x > 0].mean() if (x > 0).any() else 0:.2f}")

Print that table next to the plot. Two clusters can produce identical violin outlines while differing between 20% and 70% detection, because detection rate and expression level among expressing cells trade off against each other. If the fraction-detected distinction is the point you are making, use sc.pl.dotplot, where dot size encodes the fraction expressing and color encodes mean expression among all cells. Dot plots are strictly more informative than violins for marker panels, and we reach for them first. Violins earn their place when the shape itself matters: bimodality, a long tail, or a subpopulation hiding inside a cluster.

6. Add the statistics separately, and get them from a model

Inference should come from a model fit to the appropriate data, not from the values you happened to plot. Do not compute a t-test on the plotted values and annotate a p-value onto the figure. For bulk data, get the p-value from DESeq2 on raw counts:

dds <- DESeq(dds)
res <- results(dds, contrast = c("condition", "treated", "control"),
               alpha = 0.05)
res["HMOX1", c("log2FoldChange", "lfcSE", "pvalue", "padj")]

Then write the adjusted p-value into the caption. The violin shows the data and the model supplies the inference, and keeping them as separate objects is what makes each interpretable. For single-cell data, sc.tl.rank_genes_groups(adata, "leiden", method="wilcoxon", tie_correct=True) run on log-normalized values gives a defensible ranking, with one caveat: per-cell tests treat cells as independent replicates, so p-values are inflated when all cells come from one donor. If you have multiple donors, aggregate to pseudobulk per donor per cluster and run DESeq2 on that. This is the design multi-omics papers use when a claim needs to survive replication³⁴.

7. Cross-check the violin against an orthogonal measurement

A violin plot is one view of one assay, which limits how far it can carry a conclusion. Transcript abundance and protein abundance correlate loosely, typically in the 0.4-0.6 range across genes, so a clean RNA violin does not establish that the protein follows. Published multi-omics work treats the RNA violin as a hypothesis and confirms it with an independent layer such as immunohistochemistry, proteomics, or spatial data. The HMGB2 hepatocellular carcinoma analysis pairs single-cell expression with bulk cohorts and protein-level validation before claiming a role in the tumor microenvironment⁵, and spatial studies of endometriomas combine single-cell violins with spatially resolved metabolite and transcript measurements to confirm that a signal sits where the clustering says it does⁴.

On a personal dataset, your orthogonal layers are proteomics and blood biomarkers. If an RNA violin shows something interesting, look for the corresponding protein before you believe it. Nothing here is a clinical finding; if a result looks like it bears on your health, take the raw values and the caveats to a physician who can order a validated assay.

Common problems

Most violin plots fail in a handful of recognizable ways. Each of the following describes the symptom, its cause, and the fix.

The y-axis has no units. Label it with something specific: log2 CPM (prior.count=2), VST, or log1p(CP10K). A violin plot with a bare “Expression” axis cannot be interpreted by anyone, including you in three months.

The violin extends below zero on a log scale. You left cut at its default. Set cut=0 in seaborn and scanpy, or trim = TRUE in Seurat::VlnPlot.

Two lobes appear and disappear as you change bandwidth. They are not real. Confirm bimodality with a method that does not involve smoothing. Options include a histogram with a fixed bin count, an ECDF, or Hartigan’s dip test (diptest::dip.test in R). If the dip test on the raw values gives p > 0.1, do not describe the gene as bimodal.

Every cluster looks the same width. You combined scale="width" with standard_scale="var" and normalized away the contrast you were trying to show. Pick one normalization.

Scaled values got plotted. In scanpy, pass use_raw=True or layer="lognorm". In Seurat, use VlnPlot(obj, features = "GENE", slot = "data") rather than slot = "scale.data". A violin centered at zero with negative values is the tell.

Group sizes are wildly different and unlabeled. A cluster of 31 cells and a cluster of 4,200 cells drawn at equal width invites a false comparison. Put n= in the tick labels, every time.

Batch structure inside a group. If one condition was sequenced on a different run, a bimodal violin may be showing you two batches. Color the stripplot points by batch and look. If the lobes separate by batch, fix the design by adding batch to the DESeq2 formula or integrating with harmonypy or scvi-tools before interpreting shape.

You are plotting a gene that barely passed filtering. Check the raw counts. A gene with a maximum of 8 counts in any sample will produce a perfectly respectable-looking VST violin that represents nothing but Poisson noise. Print the raw count range in the caption for any gene below about 50 counts at its maximum.

Woolf Software builds longitudinal molecular profiles of individuals: whole-genome sequencing, RNA sequencing, proteomics, blood biomarkers, and continuous glucose data, integrated into one model of you. Build your profile.

Footnotes

1. Wenwen Shao, Zhiheng Lin, Zhikai Xiahou, et al. Single-cell RNA sequencing reveals that MYBL2 in malignant epithelial cells is involved in the development and progression of ovarian cancer. Frontiers in Immunology, 2024. https://doi.org/10.3389/fimmu.2024.1438198 ↩

2. Wen-jin Chen, Ke-qin Dong, Xiu-wu Pan, et al. Single-cell RNA-seq integrated with multi-omics reveals SERPINE2 as a target for metastasis in advanced renal cell carcinoma. Cell Death & Disease, 2023. https://doi.org/10.1038/s41419-023-05566-w ↩

3. Peng Liu, Xinpei Deng, Huamao Zhou, et al. Multi-omics analyses unravel DNA damage repair-related clusters in breast cancer with experimental validation. Frontiers in Immunology, 2023. https://doi.org/10.3389/fimmu.2023.1297180 ↩

4. Yujuan Qi, Xia Chen, Sen Zheng, et al. Single-cell and spatially resolved omics reveal transcriptional and metabolic signatures of ovarian endometriomas. Nature Communications, 2025. https://doi.org/10.1038/s41467-025-66706-8 ↩ ↩²

5. Yan-zhu Chen, Zhi-shang Meng, Zuo-lin Xiang. HMGB2 drives tumor progression and shapes the immunosuppressive microenvironment in hepatocellular carcinoma: insights from multi-omics analysis. Frontiers in Immunology, 2024. https://doi.org/10.3389/fimmu.2024.1415435 ↩