Grouped Summaries

The most common question asked of a dataset is “what is the average per group?” Mean expression per gene, counts per condition, a standard deviation per batch. Both dplyr and Polars answer it the same way: split the rows into groups, then compute one or more summaries per group. dplyr calls the second step summarise; Polars calls it agg.

Every code block runs in a container in the companion repository, on the same shared counts table used across this section. The two languages agree on every number. They print the groups in a different order, dplyr sorted by the grouping key and Polars in the order the groups first appear, which is expected and harmless.

One summary per group

Mean counts per gene.

R

library(dplyr)

expr |>
  group_by(gene) |>
  summarise(mean_counts = mean(counts), .groups = "drop")

# A tibble: 6 × 2
  gene  mean_counts
  <chr>       <dbl>
1 ACTB        3881.
2 BRCA1        182.
3 EGFR         513.
4 GAPDH       4028.
5 MYC          502.
6 TP53         504.

Python

import polars as pl

expr.group_by("gene", maintain_order=True).agg(
    pl.col("counts").mean().alias("mean_counts")
)

shape: (6, 2)
┌───────┬─────────────┐
│ gene  ┆ mean_counts │
│ ---   ┆ ---         │
│ str   ┆ f64         │
╞═══════╪═════════════╡
│ TP53  ┆ 504.125     │
│ EGFR  ┆ 512.75      │
│ MYC   ┆ 501.625     │
│ BRCA1 ┆ 181.625     │
│ GAPDH ┆ 4028.25     │
│ ACTB  ┆ 3880.625    │
└───────┴─────────────┘

The tibble rounds the means for display; the values are the same. Polars keeps the groups in first-seen order, so TP53 leads.

Several summaries at once

Mean, standard deviation, and group size together. In dplyr each is a named argument to summarise; in Polars each is an expression in agg.

R

expr |>
  group_by(gene) |>
  summarise(mean_counts = mean(counts),
            sd_counts   = sd(counts),
            n           = n(),
            .groups = "drop")

# A tibble: 6 × 4
  gene  mean_counts sd_counts     n
  <chr>       <dbl>     <dbl> <int>
1 ACTB        3881.     476.      8
2 BRCA1        182.      20.6     8
3 EGFR         513.     241.      8
4 GAPDH       4028.     262.      8
5 MYC          502.     291.      8
6 TP53         504.     101.      8

Python

expr.group_by("gene", maintain_order=True).agg(
    pl.col("counts").mean().alias("mean_counts"),
    pl.col("counts").std().alias("sd_counts"),
    pl.len().alias("n"),
)

shape: (6, 4)
┌───────┬─────────────┬────────────┬─────┐
│ gene  ┆ mean_counts ┆ sd_counts  ┆ n   │
│ ---   ┆ ---         ┆ ---        ┆ --- │
│ str   ┆ f64         ┆ f64        ┆ u32 │
╞═══════╪═════════════╪════════════╪═════╡
│ TP53  ┆ 504.125     ┆ 101.426169 ┆ 8   │
│ EGFR  ┆ 512.75      ┆ 240.643631 ┆ 8   │
│ MYC   ┆ 501.625     ┆ 290.985978 ┆ 8   │
│ BRCA1 ┆ 181.625     ┆ 20.570002  ┆ 8   │
│ GAPDH ┆ 4028.25     ┆ 261.944787 ┆ 8   │
│ ACTB  ┆ 3880.625    ┆ 476.141017 ┆ 8   │
└───────┴─────────────┴────────────┴─────┘

Both use the sample standard deviation (divides by n - 1), so sd_counts matches to full precision.

Group by more than one key

Grouping by gene and condition gives the treated versus control mean for each gene, which is the shape of a differential-expression summary.

R

expr |>
  group_by(gene, condition) |>
  summarise(mean_counts = mean(counts), .groups = "drop") |>
  arrange(gene, condition)

# A tibble: 12 × 3
   gene  condition mean_counts
   <chr> <chr>           <dbl>
 1 ACTB  control         3968.
 2 ACTB  treated         3794.
 3 BRCA1 control          171.
 4 BRCA1 treated          192
 5 EGFR  control          296
 6 EGFR  treated          730.
 7 GAPDH control         4005.
 8 GAPDH treated         4052.
 9 MYC   control          246
10 MYC   treated          757.
11 TP53  control          411.
12 TP53  treated          598.

Python

(
    expr.group_by(["gene", "condition"], maintain_order=True)
    .agg(pl.col("counts").mean().alias("mean_counts"))
    .sort(["gene", "condition"])
)

shape: (12, 3)
┌───────┬───────────┬─────────────┐
│ gene  ┆ condition ┆ mean_counts │
│ ---   ┆ ---       ┆ ---         │
│ str   ┆ str       ┆ f64         │
╞═══════╪═══════════╪═════════════╡
│ ACTB  ┆ control   ┆ 3967.5      │
│ ACTB  ┆ treated   ┆ 3793.75     │
│ BRCA1 ┆ control   ┆ 171.25      │
│ BRCA1 ┆ treated   ┆ 192.0       │
│ EGFR  ┆ control   ┆ 296.0       │
│ …     ┆ …         ┆ …           │
│ GAPDH ┆ treated   ┆ 4051.75     │
│ MYC   ┆ control   ┆ 246.0       │
│ MYC   ┆ treated   ┆ 757.25      │
│ TP53  ┆ control   ┆ 410.75      │
│ TP53  ┆ treated   ┆ 597.5       │
└───────┴───────────┴─────────────┘

EGFR jumps from a control mean near 296 to a treated mean near 730, and MYC from 246 to 757. Those are the growth genes the fixture was built to move.

Count rows per group

A grouped count needs no value column, just the group sizes.

R

count(expr, condition)

  condition  n
1   control 24
2   treated 24

Python

expr.group_by("condition", maintain_order=True).agg(pl.len().alias("n"))

shape: (2, 2)
┌───────────┬─────┐
│ condition ┆ n   │
│ ---       ┆ --- │
│ str       ┆ u32 │
╞═══════════╪═════╡
│ control   ┆ 24  │
│ treated   ┆ 24  │
└───────────┴─────┘

Twenty-four measurements per condition: four samples times six genes. Next, Joins attach gene metadata to these summaries.