R shines in data analysis thanks to the tidyverse for intuitive data wrangling and ggplot2 for high-quality visualization though very large datasets often require tools like DelayedArray or DuckDB. Python, meanwhile, leads in deep learning and scales naturally for pipelines and application development, but its data manipulation can feel less fluid and its plotting options less cohesive. In today’s post, we’ll explore how packages like plotnine, siuba, and pyjanitor help close the gap between the two ecosystems.

As described in my last post, Use uv for python and forget about anything else!, I used uv to create an environment and install all python packages on my mac laptop in less than a second. Just create a new folder for the project, open a terminal and run

1uv init
2uv add pandas siuba pyjanitor plotnine

This will create a new Python virtual environment and install the required packages.

While writing the post I needed additional packages, it was a easy as opening a terminal on the folder of my project and running uv add again.

1uv add ipykernel
2uv add setuptools
3uv add scanpy
4uv add certifi
5uv add itertools

Now in Jupyter notebook or wherever you like coding Python just load the packages as usual

1import pandas as pd
2from siuba import _, filter, mutate, summarize, group_by, arrange
3import janitor
4from plotnine import *

We will start building a ggplot with the classical R dataset mtcars. Load dataset as follows.

1from siuba.data import mtcars

/Users/alfonsosaeravila/Documents/plotnine test/.venv/lib/python3.13/site-packages/siuba/data/__init__.py:21: UserWarning: pkg_resources is deprecated as an API. See https://setuptools.pypa.io/en/latest/pkg_resources.html. The pkg_resources package is slated for removal as early as 2025-11-30. Refrain from using this package or pin to Setuptools<81.

Some data manipulation to filter and add new variables.

 1df = (
 2    mtcars
 3    >> filter(_.mpg > 20)                            # Equivalent to filter(mpg > 20)
 4    >> mutate(power_to_weight = _.hp / _.wt)         # Equivalent to mutate(power_to_weight = hp / wt)
 5    >> group_by(_.cyl)
 6    >> summarize(avg_mpg = _.mpg.mean(),
 7                 avg_power_to_weight = _.power_to_weight.mean())
 8)
 9
10print(df)

   cyl    avg_mpg  avg_power_to_weight
0    4  26.663636            37.925332
1    6  21.133333            38.153407


/Users/alfonsosaeravila/Documents/plotnine test/.venv/lib/python3.13/site-packages/siuba/dply/verbs.py:581: FutureWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.

ggplot needs the data in long format. Let's prep the data for the plot as follows.

1df_long = (
2    df
3    .pivot_longer(
4        index="cyl",
5        names_to="variable",
6        values_to="value"
7    )
8)
9print(df_long)

   cyl             variable      value
0    4              avg_mpg  26.663636
1    6              avg_mpg  21.133333
2    4  avg_power_to_weight  37.925332
3    6  avg_power_to_weight  38.153407

Now just make the plot using ggplot, you can see how similar the syntax is to R

 1(
 2    ggplot(df_long, aes(x="factor(cyl)", y="value", fill="variable"))
 3    + geom_col(position="dodge")
 4    + labs(
 5        title="Fuel efficiency and power-to-weight ratio by cylinder count",
 6        x="Cylinders",
 7        y="Average value",
 8        fill="Metric"
 9    )
10    + theme_bw()
11)

Matplotlib is building the font cache; this may take a moment.

And that is it! We have created our first ggplot plot in Python.

Now let's apply it to a single cell dataset. We will do the test creating a barplot with the mean expression of several genes by cluster, and present each gene in a single plot (ggplot facet style).

I will use the pbmc3k dataset that comes alreay pre-processed with scanpy.

1import scanpy as sc
2adata = sc.datasets.pbmc3k_processed()

WARNING: Failed to open the url with default certificates. Trying to use certifi.


100%|██████████| 23.5M/23.5M [00:04<00:00, 5.70MB/s]

Processed 3k PBMCs from 10x Genomics.
Processed using the basic tutorial Preprocessing and clustering 3k PBMCs (legacy workflow).
For preprocessing, cells are filtered out that have few gene counts or too high a percent_mito. The counts are logarithmized and only genes marked by highly_variable_genes() are retained. The obs variables n_counts and percent_mito are corrected for using regress_out(), and values are scaled and clipped by scale(). Finally, pca() and neighbors() are calculated.

We can explore the adata object.

1adata

AnnData object with n_obs × n_vars = 2638 × 1838
    obs: 'n_genes', 'percent_mito', 'n_counts', 'louvain'
    var: 'n_cells'
    uns: 'draw_graph', 'louvain', 'louvain_colors', 'neighbors', 'pca', 'rank_genes_groups'
    obsm: 'X_pca', 'X_tsne', 'X_umap', 'X_draw_graph_fr'
    varm: 'PCs'
    obsp: 'distances', 'connectivities'

and the gene expression.

1adata.to_df()

I am not sure why some expression values are negative. In any case, it will not impact the test so I just move on.

We define the top expressed genes as follows.

1import numpy as np
2scores = adata.X.sum(axis=0)
3indexes = np.argsort(-scores)[:6]
4genes = [adata.var_names[i] for i in indexes]

and extract gene expression values

1expr_df = adata[:, genes].to_df()        # cells x selected genes
2expr_df["louvain"] = adata.obs["louvain"].astype(str).values
3expr_df["cell_id"] = expr_df.index

Now we will calculate the mean expression of each gene for each louvain cluster as follows:

 1plot_df = (
 2    expr_df
 3    .pivot_longer(
 4        index=["cell_id","louvain"],
 5        names_to="gene",
 6        values_to="expr"
 7    )
 8    >> group_by(_.louvain, _.index)
 9    >> summarize(mean_expr = _.expr.mean())
10    >> arrange(_.louvain, _.index)
11)

/Users/alfonsosaeravila/Documents/plotnine test/.venv/lib/python3.13/site-packages/siuba/dply/verbs.py:581: FutureWarning: DataFrameGroupBy.apply operated on the grouping columns. This behavior is deprecated, and in a future version of pandas the grouping columns will be excluded from the operation. Either pass `include_groups=False` to exclude the groupings or explicitly select the grouping columns after groupby to silence this warning.

The last step is plotting the data.

 1(
 2    ggplot(plot_df, aes(x="index", y="mean_expr", fill = 'index'))
 3    + geom_col()
 4    + facet_wrap("~ louvain")
 5    + theme(axis_text_x=element_text(rotation=90, ha="center"), legend_position='none')
 6    + labs(
 7        x="Gene",
 8        y="Mean expression",
 9        title="Mean expression per gene, faceted by Louvain cluster"
10    )
11)

Take home messages

• Python can mimic the tidyverse workflow using siuba for dplyr-style data manipulation, pyjanitor for clean and expressive pipelines, and plotnine for ggplot2-like graphics.
• The transition from R to Python feels natural when using these packages: you keep the grammar, readability, and composability of R’s ecosystem while working inside Python’s scalable, production-oriented environment.
• These tools scale to real bioinformatics datasets, allowing tidyverse-style wrangling and ggplot-style visualization directly on objects like AnnData—making Python a stronger end-to-end option for single-cell workflows.

Further reading

• A grammar of graphics for Python Plotnine package
• Pandas package
• Python port of dplyr and other R libraries siuba
• Clean user-friendly API pandas extension for powerful and readable data-cleaning functions pyjanitor
• R ggplot2 package
• R Tidvverse package