From raw data to final figure: a complete R workflow for bench biologists

Table of contents

The payoff
What you'll learn
Setup
Step 1: Load
Step 2: Clean
Step 3: Summarize
Step 4: Plot
Step 5: Annotate significance
Step 6: Export
The complete script
Applying it to ELISA data
What's next
Resources

This is Arc 1, Part 6 of the R for Biologists series.

Updated November 2026: statistics and error bars now use biological replicates (technical replicates averaged first), every output was regenerated by running the code, and the figures are now shown in the post.

The payoff

You've built each piece separately.

In Post 2, you loaded a qPCR CSV and explored it with head() and glimpse(). In Post 3, you cleaned it with dplyr — filtered out undetermined wells, averaged technical replicates, calculated delta Ct. In Post 4, you turned those summaries into a dot plot with error bars. In Post 5, you ran ANOVA and Tukey HSD and got your p-values.

Now we put it all together.

This post runs the complete workflow — raw CSV to publication-ready figure with significance brackets — in one script. By the end you'll have a template you can adapt for any lab dataset.

What you'll learn

By the end of this post, you'll be able to:

Run a complete analysis from raw CSV to exported figure in a single R script
Add significance brackets and p-value annotations to a ggplot2 figure with ggsignif
Export at 300 dpi for journal submission
Apply the same 6-step workflow to a new dataset (ELISA)

Setup

Install ggsignif — the one new package in this post. Everything else you already have from the earlier posts:

install.packages("ggsignif")  # new in this post — skip if already installed

Load all four packages:

library(dplyr)
library(readr)
library(ggplot2)
library(ggsignif)

Download qpcr_long.csv if you don't have it — same file from Post 2 onward.

Step 1: Load

data <- read_csv("qpcr_long.csv")
glimpse(data)

Rows: 90
Columns: 8
$ sample_id      <chr> "S01_T1", "S01_T2", "S02_T1", "S02_T2", "S03_T1", "S03_T2", "S04_…
$ group          <chr> "Control", "Control", "Control", "Control", "Control", "Control",…
$ biological_rep <dbl> 1, 1, 2, 2, 3, 3, 1, 1, 2, 2, 3, 3, 1, 1, 2, 2, 3, 3, 1, 1, 2, 2,…
$ technical_rep  <dbl> 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2, 1, 2,…
$ gene           <chr> "GAPDH", "GAPDH", "GAPDH", "GAPDH", "GAPDH", "GAPDH", "GAPDH", "G…
$ ct_value       <dbl> 21.15, 21.22, 20.87, 20.94, 21.38, 21.31, 21.09, 21.14, 21.31, 21…
$ plate_id       <chr> "P01", "P01", "P01", "P01", "P01", "P01", "P01", "P01", "P02", "P…
$ date           <date> 2026-01-08, 2026-01-08, 2026-01-08, 2026-01-08, 2026-01-08, 2026…

Same 90-row file from Post 2. One row per measurement — one gene, one sample, one technical replicate.

Step 2: Clean

Drop undetermined wells (Ct ≥ 35) and the ACTB reference gene — we'll use GAPDH only. Then keep the columns we actually need:

clean <- data |>
  filter(ct_value < 35) |>
  filter(!gene %in% c("ACTB")) |>
  select(sample_id, group, biological_rep, gene, ct_value)

nrow(clean)

[1] 72

72 rows: 4 genes × 3 groups × 3 biological replicates × 2 technical replicates. All Ct values in this dataset fall below 35, so the first filter keeps everything; its value is as a safeguard for real data that might have undetermined wells.

Step 3: Summarize

The same idea as Post 3, in two stages. First average the two technical replicates of each sample, so you have one value per biological replicate. Then summarize across the biological replicates:

# Stage 1: one value per biological replicate
bio_means <- clean |>
  group_by(group, gene, biological_rep) |>
  summarise(ct_value = mean(ct_value), .groups = "drop")

# Stage 2: mean and SD across biological replicates
summary_stats <- bio_means |>
  group_by(group, gene) |>
  summarise(
    mean_ct = mean(ct_value),
    sd_ct   = sd(ct_value),
    n       = n(),
    .groups = "drop"
  )

print(summary_stats)

# A tibble: 12 × 5
   group    gene  mean_ct sd_ct     n
   <chr>    <chr>   <dbl> <dbl> <int>
 1 Control  GAPDH    21.1 0.223     3
 2 Control  IL10     32.0 0.340     3
 3 Control  IL6      30.9 0.418     3
 4 Control  TNF      29.8 0.265     3
 5 LPS_10ng GAPDH    21.2 0.235     3
 6 LPS_10ng IL10     27.2 0.380     3
 7 LPS_10ng IL6      23.9 0.342     3
 8 LPS_10ng TNF      24.3 0.379     3
 9 LPS_1ng  GAPDH    21.1 0.224     3
10 LPS_1ng  IL10     29.4 0.478     3
11 LPS_1ng  IL6      27.3 0.420     3
12 LPS_1ng  TNF      27.2 0.203     3

Notice n = 3, not 6. The two wells of a sample are the same cDNA measured twice, so they are not independent observations. Your error bars and your statistics should reflect the 3 biological replicates you actually have. (Why this matters is covered in Post 5.)

We'll focus on IL-6 for the figure. The dose-response is already visible in the numbers: Control 30.9 → LPS_1ng 27.3 → LPS_10ng 23.9. Each ~3.5 Ct drop is roughly an 11-fold increase in expression (2^3.5).

Step 4: Plot

Filter to IL-6 and set the group order as a factor. We make two tables: the summary (for the mean and error bar) and the per-replicate values (so readers can see the individual samples). This is the scaffold from Post 4, plus one layer for the replicates:

grp <- c("Control", "LPS_1ng", "LPS_10ng")

il6_summary <- summary_stats |>
  filter(gene == "IL6") |>
  mutate(group = factor(group, levels = grp))

il6_bio <- bio_means |>
  filter(gene == "IL6") |>
  mutate(group = factor(group, levels = grp))

p <- ggplot(il6_summary, aes(x = group, y = mean_ct, color = group)) +
  geom_point(
    data = il6_bio, aes(y = ct_value),
    size = 2, alpha = 0.5,
    position = position_jitter(width = 0.08, seed = 1)
  ) +
  geom_point(size = 4) +
  geom_errorbar(
    aes(ymin = mean_ct - sd_ct, ymax = mean_ct + sd_ct),
    width = 0.2
  ) +
  labs(
    title = "IL-6 expression by treatment group",
    x     = NULL,
    y     = "Mean Ct (lower = higher expression)"
  ) +
  theme_classic() +
  theme(legend.position = "none")

p

Dot plot of IL-6 Ct values for Control, LPS_1ng and LPS_10ng: large points show group means with SD error bars, small faded points show the three biological replicates per group. Ct falls from about 31 in Control to about 27 at 1 ng and 24 at 10 ng LPS.

Figure 1. The plot p produces: group mean ± SD (large points) over the three biological replicates (small points). Showing the individual replicates lets readers see that n = 3.

The dose-response is visible — Control has the highest Ct (least IL-6), LPS_10ng the lowest. What's missing: significance annotations to tell readers which differences are statistically meaningful.

Step 5: Annotate significance

ggsignif adds significance brackets as a ggplot2 geom layer — the same way geom_point() or geom_errorbar() adds dots and bars. You specify which groups to compare and what label to show; it draws the bracket.

First, get the p-values from Tukey HSD (the same test from Post 5), run on the biological replicates:

anova_model <- aov(ct_value ~ group, data = il6_bio)
TukeyHSD(anova_model)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = ct_value ~ group, data = il6_bio)

$group
                      diff       lwr       upr    p adj
LPS_1ng-Control  -3.531667 -4.521332 -2.542002 8.50e-05
LPS_10ng-Control -6.968333 -7.957998 -5.978668 1.60e-06
LPS_10ng-LPS_1ng -3.436667 -4.426332 -2.447002 9.94e-05

All three comparisons are p < 0.001 — we'll use *** for all three. Now add the brackets:

p_annotated <- p +
  geom_signif(
    comparisons = list(
      c("Control",  "LPS_1ng"),
      c("LPS_1ng",  "LPS_10ng"),
      c("Control",  "LPS_10ng")
    ),
    annotations = c("***", "***", "***"),
    y_position  = c(32.0, 33.0, 34.5),
    tip_length  = 0.01,
    color       = "black"
  ) +
  expand_limits(y = 36)

p_annotated

The same IL-6 dot plot with three significance brackets above the points: Control vs LPS_1ng, LPS_1ng vs LPS_10ng, and Control vs LPS_10ng, each labelled with three asterisks.

Figure 2. The finished figure, p_annotated: each bracket marks a Tukey HSD comparison with p < 0.001.

geom_signif() takes three key arguments:

comparisons — a list of two-element vectors naming the groups to connect. Group names must match your x-axis factor levels exactly.
annotations — the label for each bracket. Use "***" for p < 0.001, "**" for p < 0.01, "*" for p < 0.05, "ns" for not significant. Or pass the numeric p-value as a string: "p = 0.03".
y_position — where to place the horizontal bar of each bracket, in data coordinates. Stack them: the innermost bracket lowest, the outermost highest.

The expand_limits(y = 36) line pushes the y axis up to make room for all three brackets. Adjust this value if brackets run off the top of your figure.

Step 6: Export

ggsave("il6_figure.png",  plot = p_annotated, dpi = 300, width = 5, height = 5)
ggsave("il6_figure.tiff", plot = p_annotated, dpi = 300, width = 5, height = 5)

300 dpi is the minimum most journals require for figures. Export both PNG (for previewing and presentations) and TIFF (for submission — many journals specify TIFF or EPS). Adjust width and height (in inches) to match your journal's column width guidelines.

The complete script

Here's the full workflow in one clean script. Save it as qpcr_analysis.R, swap in your own file path and column names, and it will run start to finish:

# === qPCR analysis: complete Arc 1 workflow ===

# Packages
library(dplyr)
library(readr)
library(ggplot2)
library(ggsignif)

# Step 1: Load
data <- read_csv("qpcr_long.csv")

# Step 2: Clean
clean <- data |>
  filter(ct_value < 35) |>
  filter(!gene %in% c("ACTB")) |>
  select(sample_id, group, biological_rep, gene, ct_value)

# Step 3: Summarize (technical reps -> biological reps -> group summary)
bio_means <- clean |>
  group_by(group, gene, biological_rep) |>
  summarise(ct_value = mean(ct_value), .groups = "drop")

summary_stats <- bio_means |>
  group_by(group, gene) |>
  summarise(
    mean_ct = mean(ct_value),
    sd_ct   = sd(ct_value),
    n       = n(),
    .groups = "drop"
  )

# Step 4: Plot — IL-6 dot plot with replicates and error bars
grp <- c("Control", "LPS_1ng", "LPS_10ng")

il6_summary <- summary_stats |>
  filter(gene == "IL6") |>
  mutate(group = factor(group, levels = grp))

il6_bio <- bio_means |>
  filter(gene == "IL6") |>
  mutate(group = factor(group, levels = grp))

p <- ggplot(il6_summary, aes(x = group, y = mean_ct, color = group)) +
  geom_point(
    data = il6_bio, aes(y = ct_value),
    size = 2, alpha = 0.5,
    position = position_jitter(width = 0.08, seed = 1)
  ) +
  geom_point(size = 4) +
  geom_errorbar(
    aes(ymin = mean_ct - sd_ct, ymax = mean_ct + sd_ct),
    width = 0.2
  ) +
  labs(
    title = "IL-6 expression by treatment group",
    x     = NULL,
    y     = "Mean Ct (lower = higher expression)"
  ) +
  theme_classic() +
  theme(legend.position = "none")

# Step 5: Significance annotations
anova_model <- aov(ct_value ~ group, data = il6_bio)
TukeyHSD(anova_model)  # check p-values before annotating

p_annotated <- p +
  geom_signif(
    comparisons = list(
      c("Control",  "LPS_1ng"),
      c("LPS_1ng",  "LPS_10ng"),
      c("Control",  "LPS_10ng")
    ),
    annotations = c("***", "***", "***"),
    y_position  = c(32.0, 33.0, 34.5),
    tip_length  = 0.01,
    color       = "black"
  ) +
  expand_limits(y = 36)

p_annotated

# Step 6: Export
ggsave("il6_figure.png",  plot = p_annotated, dpi = 300, width = 5, height = 5)
ggsave("il6_figure.tiff", plot = p_annotated, dpi = 300, width = 5, height = 5)

Applying it to ELISA data

The six steps above aren't specific to qPCR. Here's the same workflow on a plate reader dataset — OD450 absorbance readings from a cytokine ELISA, three conditions, three biological replicates each, every sample read in duplicate wells.

Download elisa_data.csv, then run steps 1–3 and the test:

# Step 1: Load
elisa <- read_csv("elisa_data.csv")

# Step 2: Clean
elisa_clean <- elisa |>
  filter(od_value > 0.05) |>
  select(sample_id, group, biological_rep, od_value) |>
  mutate(group = factor(group, levels = c("Control", "Drug_1uM", "Drug_10uM")))

# Step 3: Summarize (duplicate wells -> biological reps -> group summary)
elisa_bio <- elisa_clean |>
  group_by(group, biological_rep) |>
  summarise(od_value = mean(od_value), .groups = "drop")

elisa_summary <- elisa_bio |>
  group_by(group) |>
  summarise(
    mean_od = mean(od_value),
    sd_od   = sd(od_value),
    n       = n(),
    .groups = "drop"
  )

anova_elisa <- aov(od_value ~ group, data = elisa_bio)
TukeyHSD(anova_elisa)

  Tukey multiple comparisons of means
    95% family-wise confidence level

Fit: aov(formula = od_value ~ group, data = elisa_bio)

$group
                        diff       lwr       upr p adj
Drug_1uM-Control   0.2896667 0.2633816 0.3159517 2e-07
Drug_10uM-Control  0.6946667 0.6683816 0.7209517 0e+00
Drug_10uM-Drug_1uM 0.4050000 0.3787149 0.4312851 0e+00

Both doses differ from Control with p < 0.001 (R prints 0e+00 when the p-value is smaller than it displays). Now steps 4–6, plot, annotate and export:

ggplot(elisa_summary, aes(x = group, y = mean_od, color = group)) +
  geom_point(
    data = elisa_bio, aes(y = od_value),
    size = 2, alpha = 0.5,
    position = position_jitter(width = 0.08, seed = 1)
  ) +
  geom_point(size = 4) +
  geom_errorbar(
    aes(ymin = mean_od - sd_od, ymax = mean_od + sd_od),
    width = 0.2
  ) +
  geom_signif(
    comparisons = list(c("Control", "Drug_1uM"), c("Control", "Drug_10uM")),
    annotations = c("***", "***"),
    y_position  = c(1.0, 1.1),
    tip_length  = 0.01,
    color       = "black"
  ) +
  expand_limits(y = 1.2) +
  labs(x = NULL, y = "Absorbance (OD 450 nm)") +
  theme_classic() +
  theme(legend.position = "none")

# Step 6: Export
ggsave("elisa_figure.png", dpi = 300, width = 4, height = 5)

ELISA dot plot of OD450 absorbance for Control, Drug_1uM and Drug_10uM, rising from about 0.20 to 0.49 to 0.89, with significance brackets comparing each drug dose to Control.

Figure 3. The ELISA version of the same figure. Only the column names, filter threshold and axis label changed.

The column names changed (ct_value → od_value). The filter threshold changed (< 35 → > 0.05). The axis label changed. The rest is identical.

The steps don't change. Only the column names do.

What's next

That's Arc 1 complete.

You've gone from zero R to a full reproducible analysis pipeline: loading data, cleaning it, summarizing it, plotting it, testing it, and annotating the figure. If you've been following along with your own data, you should have a working analysis.R script by now.

Arc 2 picks up from here and moves into public datasets and bioinformatics workflows — downloading gene expression data from GEO, normalization, DESeq2, pathway enrichment. The tidyverse syntax, ggplot2 patterns, and statistical thinking from Arc 1 all carry straight over.

→ Next: How to download public gene expression datasets from GEO in R

← Previous: T-tests and ANOVA in R: lab stats without GraphPad Prism

What dataset have you been waiting to try this on? Drop it in the comments.

Resources

Resource	What it is	Link
`ggsignif`	Significance brackets for ggplot2	cran.r-project.org/package=ggsignif
`ggplot2`	Visualization (covered in Post 4)	ggplot2.tidyverse.org
`dplyr`	Data wrangling (covered in Post 3)	dplyr.tidyverse.org
`qpcr_long.csv`	qPCR dataset used throughout Arc 1	Download
`elisa_data.csv`	ELISA dataset used in the outro	Download