Computed Columns

The core verbs get you a subset of the data; the real work is the formulas you compute on it. Raw counts become log values, normalized values, and standardized scores before any model or plot. dplyr does this with mutate, Polars with with_columns. The interesting cases are the grouped ones, where the formula needs a per-sample total or a per-gene mean but you still want one row per measurement. dplyr handles that with group_by before mutate; Polars with a window expression, .over(...).

Every code block runs in a container in the companion repository. The two languages agree on every value.

A simple formula

Raw RNA-seq counts span orders of magnitude, so they are almost always log-transformed. log2(counts + 1) is the usual choice; the + 1 keeps zeros finite.

R

expr <- expr |> mutate(log2_counts = log2(counts + 1))
head(expr)

  sample_id condition  gene counts log2_counts
1       S01   control  TP53    393    8.622052
2       S01   control  EGFR    294    8.204571
3       S01   control   MYC    247    7.954196
4       S01   control BRCA1    195    7.614710
5       S01   control GAPDH   3939   11.943980
6       S01   control  ACTB   3117   11.606405

Python

expr = expr.with_columns((pl.col("counts") + 1).log(2).alias("log2_counts"))
expr.head()

shape: (5, 5)
┌───────────┬───────────┬───────┬────────┬─────────────┐
│ sample_id ┆ condition ┆ gene  ┆ counts ┆ log2_counts │
│ ---       ┆ ---       ┆ ---   ┆ ---    ┆ ---         │
│ str       ┆ str       ┆ str   ┆ i64    ┆ f64         │
╞═══════════╪═══════════╪═══════╪════════╪═════════════╡
│ S01       ┆ control   ┆ TP53  ┆ 393    ┆ 8.622052    │
│ S01       ┆ control   ┆ EGFR  ┆ 294    ┆ 8.204571    │
│ S01       ┆ control   ┆ MYC   ┆ 247    ┆ 7.954196    │
│ S01       ┆ control   ┆ BRCA1 ┆ 195    ┆ 7.61471     │
│ S01       ┆ control   ┆ GAPDH ┆ 3939   ┆ 11.94398    │
└───────────┴───────────┴───────┴────────┴─────────────┘

A grouped formula: CPM

Counts per million normalize each measurement by its sample’s total, so samples sequenced to different depths become comparable. The formula divides by a per-sample sum, which means grouping by sample_id first. In dplyr that is group_by then mutate; in Polars it is .sum().over("sample_id") inside the expression, which computes the total per sample without collapsing the rows.

R

expr <- expr |>
  group_by(sample_id) |>
  mutate(cpm = counts / sum(counts) * 1e6) |>
  ungroup()
head(expr)

# A tibble: 6 × 6
  sample_id condition gene  counts log2_counts     cpm
  <chr>     <chr>     <chr>  <int>       <dbl>   <dbl>
1 S01       control   TP53     393        8.62  48015.
2 S01       control   EGFR     294        8.20  35919.
3 S01       control   MYC      247        7.95  30177.
4 S01       control   BRCA1    195        7.61  23824.
5 S01       control   GAPDH   3939       11.9  481246.
6 S01       control   ACTB    3117       11.6  380819.

Python

expr = expr.with_columns(
    (pl.col("counts") / pl.col("counts").sum().over("sample_id") * 1e6).alias("cpm")
)
expr.head()

shape: (5, 6)
┌───────────┬───────────┬───────┬────────┬─────────────┬──────────────┐
│ sample_id ┆ condition ┆ gene  ┆ counts ┆ log2_counts ┆ cpm          │
│ ---       ┆ ---       ┆ ---   ┆ ---    ┆ ---         ┆ ---          │
│ str       ┆ str       ┆ str   ┆ i64    ┆ f64         ┆ f64          │
╞═══════════╪═══════════╪═══════╪════════╪═════════════╪══════════════╡
│ S01       ┆ control   ┆ TP53  ┆ 393    ┆ 8.622052    ┆ 48014.660965 │
│ S01       ┆ control   ┆ EGFR  ┆ 294    ┆ 8.204571    ┆ 35919.364692 │
│ S01       ┆ control   ┆ MYC   ┆ 247    ┆ 7.954196    ┆ 30177.153329 │
│ S01       ┆ control   ┆ BRCA1 ┆ 195    ┆ 7.61471     ┆ 23824.068418 │
│ S01       ┆ control   ┆ GAPDH ┆ 3939   ┆ 11.94398    ┆ 481246.18204 │
└───────────┴───────────┴───────┴────────┴─────────────┴──────────────┘

Each sample’s CPM values sum to exactly one million, which is the point of the normalization and a quick way to check you grouped by the right column.

A per-gene z-score

Standardizing each gene across samples, subtract the gene’s mean and divide by its standard deviation, puts every gene on the same scale. Same grouped pattern, now over gene. MYC, a growth gene the treatment moves, shows clearly negative z-scores in the four control samples and positive ones in the treated.

R

expr <- expr |>
  group_by(gene) |>
  mutate(z = (counts - mean(counts)) / sd(counts)) |>
  ungroup()

expr |> filter(gene == "MYC") |> select(sample_id, condition, counts, z)

# A tibble: 8 × 4
  sample_id condition counts      z
  <chr>     <chr>      <int>  <dbl>
1 S01       control      247 -0.875
2 S02       control      243 -0.889
3 S03       control      228 -0.940
4 S04       control      266 -0.810
5 S05       treated      810  1.06
6 S06       treated      531  0.101
7 S07       treated      854  1.21
8 S08       treated      834  1.14

Python

expr = expr.with_columns(
    ((pl.col("counts") - pl.col("counts").mean().over("gene"))
     / pl.col("counts").std().over("gene")).alias("z")
)

expr.filter(pl.col("gene") == "MYC").select("sample_id", "condition", "counts", "z")

shape: (8, 4)
┌───────────┬───────────┬────────┬───────────┐
│ sample_id ┆ condition ┆ counts ┆ z         │
│ ---       ┆ ---       ┆ ---    ┆ ---       │
│ str       ┆ str       ┆ i64    ┆ f64       │
╞═══════════╪═══════════╪════════╪═══════════╡
│ S01       ┆ control   ┆ 247    ┆ -0.875042 │
│ S02       ┆ control   ┆ 243    ┆ -0.888789 │
│ S03       ┆ control   ┆ 228    ┆ -0.940337 │
│ S04       ┆ control   ┆ 266    ┆ -0.809747 │
│ S05       ┆ treated   ┆ 810    ┆ 1.059759  │
│ S06       ┆ treated   ┆ 531    ┆ 0.10095   │
│ S07       ┆ treated   ┆ 854    ┆ 1.210969  │
│ S08       ┆ treated   ┆ 834    ┆ 1.142237  │
└───────────┴───────────┴────────┴───────────┘

Both use the sample standard deviation, so the z-scores match to full precision. Next, Conditional Logic turns these values into labels.