6 Session 6: Genome-wide association: thousands of linear models and one confounder

6.1 Overview

Quick poll (Poll Everywhere)

Biological question

In the same growth chamber, some wild Arabidopsis plants flower after 6 weeks, others after 4 months. Which regions of their genome explain the difference, and how do we avoid being fooled by where the plants come from?

In Session 4 you ran one linear model. In Session 5 you ran one test per molecule, 34 then 2,500 tests, and corrected the p-values. Today you run more than 20,000 linear models, one per position of the genome where the plants differ: a genome-wide association study (GWAS). The data are real and public: about 1,000 wild lines of Arabidopsis thaliana collected across Europe and Asia, sequenced by the 1001 Genomes project, and grown side by side to measure when they flower.

You will see that the naive analysis finds “significant” positions everywhere, and why: plants from the same region share both their DNA and their flowering behaviour. The PCA of Session 5, run on the DNA, finds this population structure and lets you correct for it. What remains points to a famous flowering gene. Then you do the same alone on a second trait, the sodium in the leaves, where the story is different.

How you work and submit today

As in Session 5, everything happens in WORKBOOK.qmd, and Render makes WORKBOOK.pdf, which you submit. Blocks A to E are guided; block F (leaf sodium) you do on your own. Every block counts in the grade. Submit on Moodle by Sunday evening (1 November).

New today: some chunks take 20 to 60 seconds (the PCA, and each scan of all the SNPs). Wait until the red stop sign at the top of the Console disappears before running the next one.

6.1.1 Learning outcomes

By the end of this session, you should be able to:

explain what a SNP is, and how a genotype table is coded (0/1);
test one SNP with the linear model of Session 4, then all SNPs with a function and apply();
explain why thousands of tests need a stricter threshold (Bonferroni), and draw a Manhattan plot;
read a QQ plot and the inflation factor λ, and say when a scan cannot be trusted;
explain how population structure creates false associations, find it with a PCA of the genotypes, and correct for it with principal components as covariates;
go from a peak to a list of candidate genes, and say what a GWAS can and cannot prove.

6.1.2 Session plan

Time	Block	What happens
14:15-14:20	Retrieval	Session 5 quiz.
14:20-14:35	Theory 1	Natural variation, SNPs, one linear model per SNP, Bonferroni, the Manhattan plot.
14:35-14:50	Theory 2	The confounder: population structure, the QQ plot, the PCA correction; from a peak to a gene.
14:50-15:00	Practical briefing	What today’s workbook will ask of you.
15:00-15:15	Practical A	Same project habits, new data.
15:15-15:35	Practical B	The accessions: a map, flowering time per group, without statistics.
15:35-15:45	Break
15:45-16:10	Practical C	One SNP at a time: a SNP that does nothing, and a SNP that lies.
16:10-16:35	Practical D	Guided: one linear model per SNP, the Manhattan and QQ plots.
16:35-17:10	Practical E	Guided: the population structure, the corrected GWAS, the gene.
17:10-17:15	Wrap-up	Render, and how block F (leaf sodium) works.

6.2 Theory (≈45 min)

6.2.1 1. The question: natural variation

Individuals of one species differ, and part of these differences is written in their DNA. A GWAS asks: at which positions of the genome do the individuals with a high value of a trait tend to carry a different letter than those with a low value?

Arabidopsis thaliana is ideal for this. It fertilises itself, so a wild line, an accession, is almost completely homozygous: its seeds are genetically identical to the parent plant. An accession can be sequenced once and grown as many times as needed, in as many conditions as needed. The 1001 Genomes project (2016) sequenced 1,135 accessions and found about 10 million positions where they differ.

6.2.2 2. The data: SNPs coded 0 and 1

A SNP (single nucleotide polymorphism) is one position of the genome where some accessions carry one letter and others another, for example A or G. Today’s genotype table keeps more than 20,000 SNPs, about one every 5,000 bases, coded:

0: the accession has the same letter as the reference genome (the lab accession Col-0);
1: it has the other letter.

0 does not mean “normal” and 1 does not mean “mutant”: 0 only means “like Col-0”. Each accession is homozygous, so there are only two values (in humans, who carry two different copies, the coding is 0, 1 or 2).

Only SNPs whose rarer letter is carried by at least 5% of the accessions were kept (the minor allele frequency, MAF ≥ 5%): a letter carried by 3 plants out of 1,000 would be tested on 3 plants only, and one odd plant would be enough to make it look important.

The table has one row per SNP, as in the VCF format used in genomics, and one column per accession. In R, you turn it into a matrix of 0 and 1, as you did with the metabolites in Session 5.

6.2.3 3. One linear model per SNP

For one SNP, the question is the one of Session 4: does the flowering time differ between the plants with allele 0 and those with allele 1?

lm(flowering ~ allele)

The Estimate of allele is the difference between the two means, in days; its p-value says whether the difference is larger than chance. A GWAS repeats this model for every SNP and keeps the p-values: thousands of models, each a fraction of a second.

Many tests. At the 5% level, 1 test in 20 is “significant” by chance alone: with thousands of SNPs, that makes hundreds of false positives, even if no SNP did anything. GWAS use the Bonferroni threshold: 5% divided by the number of tests. Only SNPs below it are called significant. It is stricter than the FDR correction of Session 5 (it controls the risk of even one false positive in the whole genome), and it is the usual choice in GWAS.

The Manhattan plot shows all the tests at once: the position of each SNP along the five chromosomes on the x axis, −log10(p) on the y axis (p = 0.01 becomes 2, 10⁻⁶ becomes 6: high means strong evidence), and a line at the Bonferroni threshold. A real association looks like a tower: not one point, but several neighbouring SNPs raised together, because neighbouring SNPs are inherited together (linkage disequilibrium).

6.2.4 4. Can the scan be trusted? The QQ plot

If no SNP had any effect, the p-values would be spread evenly between 0 and 1. The QQ plot compares the observed p-values (sorted) with these expected ones, both on the −log10 scale:

points on the diagonal, with a few rising at the top right only: a healthy scan, with a few real associations;
the whole curve above the diagonal, almost from the start: too many small p-values everywhere, a sign of a hidden problem.

The inflation factor λ sums up the QQ plot in one number: λ = 1 means no inflation; above about 1.1, be suspicious. It is computed from the median SNP, so a few real peaks do not change it.

6.2.5 5. The confounder: population structure

Plants from the same region share ancestry: over most of the genome, two Swedish accessions are more alike than a Swedish and a Spanish one. Thousands of letters are therefore frequent in Sweden and rare elsewhere, for purely historical reasons. If the trait also differs between regions (say, if the plants of one region flower late), then every SNP with a letter typical of that region looks associated with late flowering. The SNP did nothing: it only tells where the plant comes from. Ancestry is a confounder, like a place on the bench could have been one in the soybean experiment of Session 4, if the plants had not been randomised. In nature, nothing is randomised.

The fix. The PCA of Session 5, run on the genotype matrix (accessions in rows, SNPs in columns), finds the main axes of genetic differences between accessions: they follow ancestry, and geography. Adding the first principal components to every model,

lm(flowering ~ allele + pcs)

changes the question to: among accessions of similar ancestry, do those with allele 1 flower later? This is what adding heat next to water did in Session 4: an effect “adjusted for” another variable. Professionals go one step further with mixed models and a kinship matrix (the genetic similarity of every pair of plants); we stay with the PCs, which you already know.

6.2.6 6. From a peak to a gene

The best SNP of a tower is a marker, not necessarily the cause: it travels with the real cause, which may be another SNP, a deletion, or anything nearby. In Arabidopsis, this travelling together usually fades within about 10,000 bases, so you list the genes within ±20 kb of the best SNP, read what they do, and choose the most plausible one. Databases help: TAIR for Arabidopsis genes, UniProt for proteins (Session 3). A GWAS gives candidates; only experiments (mutants, swapping alleles between plants) prove that a gene causes the difference.

Which plot, for which question?

Your question	Plot or number	In R
Does one SNP change the trait, and by how much?	boxplot by allele, linear model	`geom_boxplot()`, `lm()`
Where in the genome are the associations?	Manhattan plot	`geom_point()` + `facet_grid()`
Can I trust the scan as a whole?	QQ plot, λ	`ppoints()`, `qchisq()`
Do the plants differ by ancestry?	PCA of the genotypes	`prcomp()`
Which genes sit under a peak?	the gene table, ±20 kb	`filter()`

6.2.7 7. Worked example: the 1001 Genomes plants

Flowering time (1001 Genomes Consortium, 2016): each accession was grown in growth chambers at a constant 10 °C or 16 °C, 16 hours of light, without any winter, 4 plants per accession. The trait is the number of days from sowing to the first open flower (the mean of the 4 plants).
Leaf sodium (Baxter et al., 2010): plants of about 150 of these accessions were grown in the same soil, and the sodium in their leaves measured by mass spectrometry. Too much sodium is toxic for plants: it is the main problem of salty soils.
Genetic groups: the 1001 Genomes project grouped the accessions by ancestry (north_sweden, spain, relict, …); accessions with mixed ancestry are admixed. You will use these groups only to look at the results, never in the GWAS itself.

6.2.8 Practical briefing

Blocks A-E in pairs, block F on your own. Block B meets the plants without any statistics: where they come from, and when they flower. Block C tests two SNPs by hand: one that does nothing, and one that looks important but lies. Block D tests all the SNPs, flowering at 16 °C; block E finds the structure with a PCA, corrects the GWAS, and looks up the gene under the best peak. Block F does the same alone for leaf sodium: same code, a different answer. Your submission is the whole rendered workbook: every block counts.

6.3 Your R toolbox

New tools today, on top of the Session 4 and Session 5 toolboxes (still valid: see the Session 4 R toolbox).

A function packs a few lines of code under a name, so that you can run them again on other data:

snp_p_value <- function(allele) {
  fit <- lm(flowering ~ allele)
  summary(fit)$coefficients["allele", "Pr(>|t|)"]
}

allele is the argument: the data you give the function each time you call it;
the last line is what the function returns: here the p-value of the row allele in the table of coefficients of summary(fit);
flowering is not an argument: the function finds it in your session, so create it first.

apply() runs a function on every row (or column) of a matrix:

apply(genotype_matrix, 1, snp_p_value)   # 1 = each row, 2 = each column

It returns one value per row: here one p-value per SNP, in the order of the rows.

Common mistake. object 'flowering' not found: run the line that creates flowering before apply().

Code	What it does
`genotype_matrix["chr1_920099", ]`	one row (one SNP), all columns (all accessions)
`genotype_matrix[, 1:5]`	all rows, the first 5 columns
`dim(genotype_matrix)`	number of rows, then of columns
`t(genotype_matrix)`	turns the matrix around: rows become columns
`rownames()`, `colnames()`	the names of the rows and of the columns

prcomp() wants the samples (here the accessions) in rows, hence t().

ggplot(gwas, aes(x = position / 1e6, y = -log10(p_value))) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_hline(yintercept = -log10(threshold), colour = "red", linetype = "dashed") +
  facet_grid(~ chromosome, scales = "free_x", space = "free_x")

position / 1e6: positions in millions of bases (Mb);
facet_grid(~ chromosome, ...): one panel per chromosome, side by side; scales = "free_x" gives each panel its own x axis, space = "free_x" makes each panel as wide as its chromosome is long.

qq <- gwas %>%
  arrange(p_value) %>%
  mutate(expected = -log10(ppoints(n())),
         observed = -log10(p_value))

ppoints(n()) gives n values spread evenly between 0 and 1: the p-values you would expect if no SNP had any effect;
geom_abline(intercept = 0, slope = 1) draws the diagonal.

The inflation factor:

median(qchisq(p_values, df = 1, lower.tail = FALSE)) / qchisq(0.5, df = 1)

qchisq() turns each p-value back into the test statistic it came from; λ compares the median statistic with the one expected under “no effect”.

world <- map_data("world")
ggplot() +
  geom_polygon(data = world, aes(x = long, y = lat, group = group),
               fill = "grey90", colour = "grey60") +
  geom_point(data = accessions, aes(x = longitude, y = latitude)) +
  coord_quickmap(xlim = c(-12, 90), ylim = c(30, 68))

map_data("world") gives the outlines of the countries (from the package maps);
group = group draws each country as its own shape;
coord_quickmap() keeps the proportions of a map, and xlim/ylim zoom in.

6.4 Practical (≈1 h 50)

Tested on LEGcompute

Everything used today is already installed in RStudio Server’s R (tidyverse, maps). Rendering the complete workbook takes about 1.5 to 2 minutes, because of the PCA and the four scans of all the SNPs.

6.4.1 A. Same habits, new data (15:00-15:15)

A1. The project. In the Terminal tab:

mkdir -p ~/bioinfo_course/session_06/data ~/bioinfo_course/session_06/results
cp /home/cmaslard/bioinfo_tools_course/session_06/arabidopsis_accessions.csv ~/bioinfo_course/session_06/data/
cp /home/cmaslard/bioinfo_tools_course/session_06/arabidopsis_genotypes.csv.gz ~/bioinfo_course/session_06/data/
cp /home/cmaslard/bioinfo_tools_course/session_06/arabidopsis_genes.csv ~/bioinfo_course/session_06/data/
cp /home/cmaslard/bioinfo_tools_course/session_06/WORKBOOK.qmd ~/bioinfo_course/session_06/
ls -R ~/bioinfo_course/session_06

The genotype file is compressed (.gz). Have a look at its first lines without opening it, still in the Terminal:

zcat ~/bioinfo_course/session_06/data/arabidopsis_genotypes.csv.gz | head -n 4 | cut -c 1-80

Then File → New Project → Existing Directory, pick ~/bioinfo_course/session_06, open WORKBOOK.qmd, write your name at the top and click Render once.

Checkpoint

WORKBOOK.pdf appears (RStudio opens it), and the first chunk lists the columns accession_id, name, country, latitude, longitude, genetic_group, flowering_10C, flowering_16C, leaf_sodium. The chunks further down show Error: unexpected input: those are the ___ you will fill. The zcat line prints snp_id,chromosome,position,acc_88,... then lines of 0 and 1.

read_csv() reads the compressed file directly: no need to decompress it.

6.4.2 B. The accessions, without statistics (15:15-15:35, pairs)

B1. How many, and which?

# how many accessions per genetic group?
accessions %>% count(___, sort = TRUE)

# how many accessions have a value for each trait? (!is.na() is TRUE when a value is there)
accessions %>%
  summarise(flowering_10C = sum(!is.na(flowering_10C)),
            flowering_16C = sum(!is.na(flowering_16C)),
            leaf_sodium   = sum(!is.na(___)))

Checkpoint

Two small tables: the first has one row per genetic group, the second one row with three counts.

B2. Where do they come from? A map, one point per accession:

world <- map_data("world")   # the outlines of the countries

# one point per accession, coloured by genetic group
ggplot() +
  geom_polygon(data = world, aes(x = long, y = lat, group = group),
               fill = "grey90", colour = "grey60", linewidth = 0.2) +
  geom_point(data = accessions, aes(x = longitude, y = latitude, colour = ___),
             size = 1) +
  coord_quickmap(xlim = c(-12, 90), ylim = c(30, 68)) +
  labs(x = "Longitude", y = "Latitude", colour = "Genetic group")

Checkpoint

A grey map of Europe, North Africa and Asia, with coloured points and a legend of the genetic groups. A grey map with black points means the colour is missing in aes().

B3. When do they flower? The plants grew at a constant 10 °C or 16 °C, without any winter:

# flowering at 16 °C, one box per genetic group, the groups sorted by their median
ggplot(accessions, aes(x = reorder(genetic_group, flowering_16C, FUN = median, na.rm = TRUE),
                       y = ___)) +
  geom_boxplot() +
  coord_flip() +
  labs(x = NULL, y = "Days to flowering at 16 °C")

# the mean flowering time per group, at both temperatures
accessions %>%
  group_by(genetic_group) %>%
  summarise(mean_16C = mean(flowering_16C, na.rm = TRUE),
            mean_10C = mean(___)) %>%
  arrange(desc(mean_16C))

reorder(genetic_group, flowering_16C, FUN = median, na.rm = TRUE) sorts the groups by their median flowering time, so that the boxes are in order.

B4. The missing values.

# the accessions without a flowering time at 16 °C: which groups?
accessions %>%
  filter(is.na(___)) %>%
  count(genetic_group, sort = TRUE)

# did they flower early or late at 10 °C?
accessions %>%
  group_by(missing_16C = is.na(flowering_16C)) %>%
  summarise(n = n(), mean_10C = mean(flowering_10C))

6.4.3 15:35-15:45: Short break

6.4.4 C. One SNP at a time (15:45-16:10, pairs)

C1. The genotypes.

# one row per SNP (as in a VCF file): its id, chromosome, position, then one column per accession
genotypes <- read_csv("data/arabidopsis_genotypes.csv.gz")
dim(genotypes)
genotypes[1:5, 1:8]   # the first 5 SNPs, the first 5 accessions

# the 0/1 only, as a matrix: one row per SNP, one column per accession
genotype_matrix <- genotypes %>%
  select(starts_with(___)) %>%
  as.matrix()
rownames(genotype_matrix) <- genotypes$___
dim(genotype_matrix)

# are the accessions in the same order as the rows of the table `accessions`?
all(colnames(genotype_matrix) == paste0("acc_", accessions$accession_id))

Checkpoint

dim() prints two numbers twice: rows, then columns (the table has 3 more columns than the matrix). genotypes[1:5, 1:8] shows snp_id, chromosome, position and 0/1 values. The last line must print TRUE: if it prints FALSE, stop and call the instructor.

C2. A first SNP. genotype_matrix["chr1_920099", ] takes one row of the matrix: the allele of every accession at this SNP, in the order of accessions.

# a first SNP: its allele in each accession, next to the flowering time
one_snp <- accessions %>%
  select(accession_id, genetic_group, flowering_16C) %>%
  mutate(allele = genotype_matrix["chr1_920099", ])

one_snp %>% count(___)

ggplot(one_snp, aes(x = factor(allele), y = flowering_16C)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(width = 0.2, alpha = 0.3) +
  labs(x = "Allele at chr1_920099", y = "Days to flowering at 16 °C")

# the linear model of Session 4: does the allele change the flowering time?
summary(lm(flowering_16C ~ ___, data = one_snp))

C3. A second SNP, the same way, then coloured by genetic group:

# a second SNP, the same way
trap_snp <- accessions %>%
  select(accession_id, genetic_group, flowering_16C) %>%
  mutate(allele = genotype_matrix[___, ])

summary(lm(flowering_16C ~ allele, data = trap_snp))

# the same plot, the points coloured by genetic group
ggplot(trap_snp, aes(x = factor(allele), y = flowering_16C)) +
  geom_boxplot(outlier.shape = NA) +
  geom_jitter(aes(colour = ___), width = 0.2, size = 1) +
  labs(x = "Allele at chr2_13280384", y = "Days to flowering at 16 °C")

# the share of accessions with allele 1, per genetic group (the mean of 0s and 1s)
trap_snp %>%
  group_by(genetic_group) %>%
  summarise(share_allele_1 = mean(___)) %>%
  arrange(desc(share_allele_1))

mean() of 0s and 1s is the share of 1s: mean(c(0, 1, 1, 1)) is 0.75.

C4. The same model, with the genetic group.

# the same model, with the genetic group added (as heat next to water in Session 4)
summary(lm(flowering_16C ~ allele + ___, data = trap_snp))

Checkpoint

C2 to C4 each print a lm() summary with a row allele. C4 has many more rows: one per genetic group (compared with admixed, the first in alphabetical order).

6.4.5 D. Guided: one linear model per SNP (16:10-16:35, pairs)

D1. All the SNPs. The model of C2, packed into a function, then run on every row of the matrix:

# the trait, in the same order as the columns of genotype_matrix
flowering <- accessions$flowering_16C

# the p-value of one SNP: the lm() of C2, written as a small function
snp_p_value <- function(allele) {
  fit <- lm(flowering ~ allele)
  summary(fit)$coefficients["allele", "Pr(>|t|)"]
}

# check it on the SNP of C3: the same p-value as there?
snp_p_value(genotype_matrix[___, ])

# now every SNP: apply() runs the function on each row (1) of the matrix (about 20 s)
p_naive <- apply(genotype_matrix, ___, snp_p_value)

# one row per SNP: its position and its p-value
gwas_naive <- genotypes %>%
  select(snp_id, chromosome, position) %>%
  mutate(p_value = ___)

gwas_naive %>% arrange(p_value)

# how many p-values below 0.05?
sum(gwas_naive$p_value < ___)

Checkpoint

The first line prints one p-value; apply() then works for about 20 seconds; the table gwas_naive has 4 columns (snp_id, chromosome, position, p_value) and one row per SNP. object 'flowering' not found means that the first line was not run.

D2. The threshold and the Manhattan plot.

# Bonferroni: 5% divided by the number of tests
threshold <- 0.05 / ___
threshold
-log10(threshold)
sum(gwas_naive$p_value < ___)

# the Manhattan plot: one point per SNP, one panel per chromosome
ggplot(gwas_naive, aes(x = position / 1e6, y = -log10(___))) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_hline(yintercept = -log10(threshold), colour = "red", linetype = "dashed") +
  facet_grid(~ chromosome, scales = "free_x", space = "free_x") +
  labs(title = "Flowering at 16 °C, naive GWAS", x = "Position (Mb)", y = "-log10(p)")

D3. The QQ plot and λ.

# expected: the p-values if no SNP had any effect (spread evenly between 0 and 1)
qq_naive <- gwas_naive %>%
  arrange(p_value) %>%
  mutate(expected = -log10(ppoints(n())),
         observed = -log10(___))

ggplot(qq_naive, aes(x = ___, y = ___)) +
  geom_point(size = 0.5) +
  geom_abline(intercept = 0, slope = 1, colour = "red") +
  labs(title = "QQ plot, naive GWAS", x = "Expected -log10(p)", y = "Observed -log10(p)")

# the inflation factor lambda: 1 means no inflation
lambda_naive <- median(qchisq(gwas_naive$p_value, df = 1, lower.tail = FALSE)) /
  qchisq(0.5, df = 1)
lambda_naive

Checkpoint

D2: one very small number, the same on the −log10 scale, a count, and a plot with five panels (the five chromosomes) and a red dashed line. D3: a plot with a red diagonal, and one number.

6.4.6 E. Guided: the population structure, then the corrected GWAS (16:35-17:10, pairs)

E1. The PCA of the genotypes. The same prcomp() as in Session 5, on all the accessions and all the SNPs (about 30 seconds):

# the PCA of Session 5: accessions in rows (t() turns the matrix), SNPs in columns (about 30 s)
genotype_pca <- prcomp(t(___), scale. = ___)

pca_variance <- round(100 * summary(genotype_pca)$importance[2, 1:10], 1)
pca_variance

pca_scores <- as_tibble(genotype_pca$x[, 1:2]) %>%
  bind_cols(select(accessions, genetic_group, latitude, longitude))

ggplot(pca_scores, aes(x = PC1, y = PC2, colour = ___)) +
  geom_point(size = 1) +
  labs(x = paste0("PC1 (", pca_variance[1], "%)"),
       y = paste0("PC2 (", pca_variance[___], "%)"),
       colour = "Genetic group")

# does PC1 follow geography?
cor(pca_scores$PC1, pca_scores$___, use = "complete.obs")
cor(pca_scores$PC1, pca_scores$longitude, use = "complete.obs")

Checkpoint

Ten percentages, a plot with one point per accession coloured by genetic group, and two correlations. On another computer the plot may be flipped left-right or up-down: the direction of an axis is arbitrary.

E2. How much of flowering follows ancestry?

# the first 10 PCs, as a matrix: one row per accession, one column per PC
pcs <- genotype_pca$x[, 1:___]

# how much of the flowering time do they explain? (R², the share of variation of Session 4)
summary(lm(flowering ~ ___))$r.squared

E3. The corrected GWAS. The function of D1, with the 10 PCs added to the model:

# the same function, with the 10 PCs in the model (as genetic_group in C4)
snp_p_value_pc <- function(allele) {
  fit <- lm(flowering ~ allele + ___)
  summary(fit)$coefficients["allele", "Pr(>|t|)"]
}

p_pc <- apply(genotype_matrix, 1, ___)   # about 30 s

gwas_pc <- genotypes %>%
  select(snp_id, chromosome, position) %>%
  mutate(p_value = ___)

sum(gwas_pc$p_value < ___)
lambda_pc <- median(qchisq(gwas_pc$p_value, df = 1, lower.tail = FALSE)) /
  qchisq(0.5, df = 1)
lambda_pc

# the Manhattan plot (as in D2)
ggplot(gwas_pc, aes(x = position / 1e6, y = -log10(___))) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_hline(yintercept = -log10(___), colour = "red", linetype = "dashed") +
  facet_grid(~ ___, scales = "free_x", space = "free_x") +
  labs(title = "Flowering at 16 °C, GWAS with 10 PCs", x = "Position (Mb)", y = "-log10(p)")

# the QQ plots of both scans, on one plot
qq_both <- bind_rows(naive = gwas_naive, corrected = gwas_pc, .id = "model") %>%
  group_by(model) %>%
  arrange(p_value, .by_group = TRUE) %>%
  mutate(expected = -log10(ppoints(n())),
         observed = -log10(p_value))

ggplot(qq_both, aes(x = expected, y = observed, colour = ___)) +
  geom_point(size = 0.5) +
  geom_abline(intercept = 0, slope = 1) +
  labs(title = "QQ plot, before and after the correction",
       x = "Expected -log10(p)", y = "Observed -log10(p)")

E4. The best peak, and its genes.

# the SNP of C3, before and after the correction
gwas_naive %>% filter(snp_id == "chr2_13280384")
gwas_pc %>% filter(snp_id == ___)

# the 10 best SNPs after the correction
gwas_pc %>% arrange(p_value) %>% head(10)

# the genes within 20 kb (20,000 bp) of the best SNP
genes <- read_csv("data/arabidopsis_genes.csv")
top_snp <- gwas_pc %>% slice_min(___, n = 1)
top_snp

genes %>%
  filter(chromosome == top_snp$___,
         end > top_snp$position - ___,
         start < top_snp$position + 20000)

# the same SNP in the naive scan: its p-value, and how many SNPs had a smaller one
naive_p_top <- gwas_naive %>% filter(snp_id == top_snp$___) %>% pull(p_value)
naive_p_top
sum(gwas_naive$p_value < naive_p_top)

E4, continued: the gene names on the plot. Nothing to fill: after your answer to E4, run this chunk. It draws the Manhattan plot of E3 again, with the best SNP in blue and the short names of the genes within 20 kb of it.

# the short names of the genes within 20 kb of the best SNP (the same genes as above)
top_genes <- genes %>%
  filter(chromosome == top_snp$chromosome,
         end > top_snp$position - 20000,
         start < top_snp$position + 20000,
         !is.na(symbol)) %>%
  mutate(symbol = str_replace_all(symbol, "%3B", ";"))   # the file writes ";" as %3B
top_label <- top_snp %>%
  mutate(genes = paste(top_genes$symbol, collapse = "\n"))   # one name per line

# the Manhattan plot of E3, with the best SNP in blue and the names of its genes
ggplot(gwas_pc, aes(x = position / 1e6, y = -log10(p_value))) +
  geom_point(size = 0.5, alpha = 0.5) +
  geom_hline(yintercept = -log10(threshold), colour = "red", linetype = "dashed") +
  geom_point(data = top_label, colour = "blue", size = 2) +
  geom_label(data = top_label, aes(label = genes), hjust = 0, vjust = 1, nudge_x = 1,
             size = 2.5) +
  facet_grid(~ chromosome, scales = "free_x", space = "free_x") +
  labs(title = "The genes near the best SNP (10 PCs)",
       x = "Position (Mb)", y = "-log10(p)")

Checkpoint

E3: a count, one number (λ), a Manhattan plot titled “GWAS with 10 PCs”, and a QQ plot with two colours. E4: tables of 4 columns, then a table of genes with gene_id, symbol, chromosome, start, end, description. A symbol shown as NA means the gene has no short name. E4, continued: the same Manhattan plot, with one blue point and a box of gene names next to it.

Then E5: summarise the flowering GWAS in a few sentences, with your own numbers.

6.4.7 F. On your own: leaf sodium

The same reasoning, on the leaf sodium of about 150 accessions. The workbook gives you the skeleton, with more ___ than in D and E: by now, most of them are the choices you made there, so reuse your own code.

F1: the distribution of sodium, how much the 3 first PCs explain, and the naive scan with its λ. lm() leaves out the accessions without a sodium value (NA) by itself. With 147 plants, 3 PCs are enough: check λ.
F2: the corrected scan and its Manhattan plot, saved with ggsave() and inserted into your report, as in Session 5, with your own caption and the label #fig-sodium.
F3: the genes under the best SNP, and the sodium of the plants by allele.
F4, comparing the two GWAS, and F5, one limitation.

6.4.8 Wrap-up (17:10-17:15)

Click Render (about 1.5 to 2 minutes for the whole workbook). Check that your PDF shows all the plots of blocks B to E. Block F follows the same logic as D and E. Submit WORKBOOK.pdf on Moodle by Sunday evening (1 November). If you see errors in your PDF, or if it does not render, submit in addition your WORKBOOK.qmd and your results/ folder (your saved plot).

How this is graded (10 points)

Criterion	Points	What it looks for
Reproducible workflow	2	The workbook renders (or the code is complete and in order), inside a clean project.
Working analysis	3	All the code of blocks A to F ran, and all your answers use the right numbers.
Required result	2	The corrected Manhattan plot of leaf sodium (F2): saved, inserted with a caption and `#fig-sodium`, cited in F4.
Biological interpretation	2	Your answers on the biology, block by block, and especially E5 and F4, with your own numbers.
Reflection	1	Your answers on the methods (the order of the accessions, the threshold, the QQ plot, λ, the correction) and your limitation (F5).

Every question counts, except the exit questions.

6.5 Checkpoint

You are done when all of the following are true:

you can say what a SNP is and what 0 and 1 mean in the genotype table;
you tested one SNP with lm(), then all of them with a function and apply();
you can explain the Bonferroni threshold, and read a Manhattan plot;
you can read a QQ plot and λ, and say when a scan cannot be trusted;
you can explain, with the SNP of C3, how population structure creates false associations;
you ran a PCA of the genotypes, and used its components to correct the GWAS;
you found the genes under a peak, and can say why a peak is not a proof;
you can explain why the structure misleads the GWAS of flowering time much more than that of leaf sodium.

6.6 Optional challenge

For your own interest, not graded.

A. How many PCs? Redo E3 with 3, then 20 PCs instead of 10. How do λ and the number of SNPs below the threshold change? Is the best SNP still the same?

B. FRIGIDA. FRI (chromosome 4, around 0.27 Mb) is one of the best-known flowering genes of Arabidopsis. Find the best SNP within 20 kb of it in the naive scan and in the corrected scan. What happened to it, and why could a correction for ancestry also hide a real gene?

C. Flowering at 10 °C. Redo D and E with flowering_10C. Is the best peak in the same place as at 16 °C?

D. A prettier Manhattan plot. Colour the points by chromosome, and label the best SNP of each tower with geom_text_repel() (package ggrepel, Session 5).

6.7 Going further

Korte A. and Farlow A. (2013). The advantages and limitations of trait analysis with GWAS: a review. Plant Methods. The best short review for plant biologists.
1001 Genomes Consortium (2016). 1,135 genomes reveal the global pattern of polymorphism in Arabidopsis thaliana. Cell. Where today’s SNPs and flowering times come from.
AraGWAS Catalog: published Arabidopsis GWAS, including flowering time, run with mixed models: compare with your peaks.
AraPheno: hundreds of other traits measured on the same accessions.

6.8 Before the next session

Session 7 moves from testing to predicting: deep learning. You will meet again the idea behind the cross-validation of Session 5: a model is only judged on data it has never seen.

--- title: "Session 6: Genome-wide association: thousands of linear models and one confounder" engine: knitr --- ## Overview [Download the lecture slides (PDF)](slides/06_lecture_slides.pdf) ::: callout-tip ## Quick poll (Poll Everywhere) Respond at [pe.app/cmaslard5](https://pe.app/cmaslard5) ::: ::: callout-note ## Biological question In the same growth chamber, some wild *Arabidopsis* plants flower after 6 weeks, others after 4 months. Which regions of their genome explain the difference, and how do we avoid being fooled by where the plants come from? ::: In Session 4 you ran **one** linear model. In Session 5 you ran one test per molecule, 34 then 2,500 tests, and corrected the p-values. Today you run **more than 20,000 linear models**, one per position of the genome where the plants differ: a **genome-wide association study** (GWAS). The data are real and public: about 1,000 wild lines of *Arabidopsis thaliana* collected across Europe and Asia, sequenced by the 1001 Genomes project, and grown side by side to measure when they flower. You will see that the naive analysis finds "significant" positions everywhere, and why: plants from the same region share both their DNA and their flowering behaviour. The PCA of Session 5, run on the DNA, finds this **population structure** and lets you correct for it. What remains points to a famous flowering gene. Then you do the same alone on a second trait, the sodium in the leaves, where the story is different. ::: callout-important ## How you work and submit today As in Session 5, everything happens in `WORKBOOK.qmd`, and **Render** makes `WORKBOOK.pdf`, which you submit. Blocks A to E are guided; block F (leaf sodium) you do on your own. Every block counts in the grade. Submit on Moodle **by Sunday evening (1 November)**. New today: some chunks take 20 to 60 seconds (the PCA, and each scan of all the SNPs). Wait until the red stop sign at the top of the Console disappears before running the next one. ::: ### Learning outcomes By the end of this session, you should be able to: - explain what a SNP is, and how a genotype table is coded (0/1); - test one SNP with the linear model of Session 4, then all SNPs with a function and `apply()`; - explain why thousands of tests need a stricter threshold (Bonferroni), and draw a Manhattan plot; - read a QQ plot and the inflation factor λ, and say when a scan cannot be trusted; - explain how population structure creates false associations, find it with a PCA of the genotypes, and correct for it with principal components as covariates; - go from a peak to a list of candidate genes, and say what a GWAS can and cannot prove. ### Session plan | Time | Block | What happens | |---:|---|---| | 14:15-14:20 | Retrieval | Session 5 quiz. | | 14:20-14:35 | Theory 1 | Natural variation, SNPs, one linear model per SNP, Bonferroni, the Manhattan plot. | | 14:35-14:50 | Theory 2 | The confounder: population structure, the QQ plot, the PCA correction; from a peak to a gene. | | 14:50-15:00 | Practical briefing | What today's workbook will ask of you. | | 15:00-15:15 | Practical A | Same project habits, new data. | | 15:15-15:35 | Practical B | The accessions: a map, flowering time per group, without statistics. | | 15:35-15:45 | Break | | | 15:45-16:10 | Practical C | One SNP at a time: a SNP that does nothing, and a SNP that lies. | | 16:10-16:35 | Practical D | Guided: one linear model per SNP, the Manhattan and QQ plots. | | 16:35-17:10 | Practical E | Guided: the population structure, the corrected GWAS, the gene. | | 17:10-17:15 | Wrap-up | Render, and how block F (leaf sodium) works. | ## Theory (≈45 min) ### 1. The question: natural variation Individuals of one species differ, and part of these differences is written in their DNA. A **GWAS** asks: *at which positions of the genome do the individuals with a high value of a trait tend to carry a different letter than those with a low value?* ***Arabidopsis thaliana*** is ideal for this. It fertilises itself, so a wild line, an **accession**, is almost completely homozygous: its seeds are genetically identical to the parent plant. An accession can be sequenced once and grown as many times as needed, in as many conditions as needed. The **1001 Genomes** project (2016) sequenced 1,135 accessions and found about 10 million positions where they differ. ### 2. The data: SNPs coded 0 and 1 A **SNP** (single nucleotide polymorphism) is one position of the genome where some accessions carry one letter and others another, for example `A` or `G`. Today's genotype table keeps more than 20,000 SNPs, about one every 5,000 bases, coded: - **0**: the accession has the same letter as the reference genome (the lab accession Col-0); - **1**: it has the other letter. 0 does not mean "normal" and 1 does not mean "mutant": 0 only means "like Col-0". Each accession is homozygous, so there are only two values (in humans, who carry two different copies, the coding is 0, 1 or 2). Only SNPs whose rarer letter is carried by at least 5% of the accessions were kept (the **minor allele frequency**, MAF ≥ 5%): a letter carried by 3 plants out of 1,000 would be tested on 3 plants only, and one odd plant would be enough to make it look important. The table has one row per SNP, as in the VCF format used in genomics, and one column per accession. In R, you turn it into a **matrix** of 0 and 1, as you did with the metabolites in Session 5. ### 3. One linear model per SNP For one SNP, the question is the one of Session 4: does the flowering time differ between the plants with allele 0 and those with allele 1? ```r lm(flowering ~ allele) ``` The `Estimate` of `allele` is the difference between the two means, in days; its p-value says whether the difference is larger than chance. A GWAS repeats this model for every SNP and keeps the p-values: thousands of models, each a fraction of a second. **Many tests.** At the 5% level, 1 test in 20 is "significant" by chance alone: with thousands of SNPs, that makes hundreds of false positives, even if no SNP did anything. GWAS use the **Bonferroni** threshold: 5% divided by the number of tests. Only SNPs below it are called significant. It is stricter than the FDR correction of Session 5 (it controls the risk of even one false positive in the whole genome), and it is the usual choice in GWAS. **The Manhattan plot** shows all the tests at once: the position of each SNP along the five chromosomes on the x axis, −log10(p) on the y axis (p = 0.01 becomes 2, 10⁻⁶ becomes 6: high means strong evidence), and a line at the Bonferroni threshold. A real association looks like a **tower**: not one point, but several neighbouring SNPs raised together, because neighbouring SNPs are inherited together (**linkage disequilibrium**). ### 4. Can the scan be trusted? The QQ plot If no SNP had any effect, the p-values would be spread evenly between 0 and 1. The **QQ plot** compares the observed p-values (sorted) with these expected ones, both on the −log10 scale: - points on the diagonal, with a few rising at the top right only: a healthy scan, with a few real associations; - the whole curve above the diagonal, almost from the start: too many small p-values **everywhere**, a sign of a hidden problem. The **inflation factor λ** sums up the QQ plot in one number: λ = 1 means no inflation; above about 1.1, be suspicious. It is computed from the median SNP, so a few real peaks do not change it. ### 5. The confounder: population structure Plants from the same region share ancestry: over most of the genome, two Swedish accessions are more alike than a Swedish and a Spanish one. Thousands of letters are therefore frequent in Sweden and rare elsewhere, for purely historical reasons. If the trait **also** differs between regions (say, if the plants of one region flower late), then every SNP with a letter typical of that region looks associated with late flowering. The SNP did nothing: it only tells where the plant comes from. Ancestry is a **confounder**, like a place on the bench could have been one in the soybean experiment of Session 4, if the plants had not been randomised. In nature, nothing is randomised. **The fix.** The PCA of Session 5, run on the genotype matrix (accessions in rows, SNPs in columns), finds the main axes of genetic differences between accessions: they follow ancestry, and geography. Adding the first principal components to every model, ```r lm(flowering ~ allele + pcs) ``` changes the question to: *among accessions of similar ancestry, do those with allele 1 flower later?* This is what adding `heat` next to `water` did in Session 4: an effect "adjusted for" another variable. Professionals go one step further with **mixed models** and a kinship matrix (the genetic similarity of every pair of plants); we stay with the PCs, which you already know. ### 6. From a peak to a gene The best SNP of a tower is a **marker**, not necessarily the cause: it travels with the real cause, which may be another SNP, a deletion, or anything nearby. In *Arabidopsis*, this travelling together usually fades within about 10,000 bases, so you list the genes within ±20 kb of the best SNP, read what they do, and choose the most plausible one. Databases help: [TAIR](https://www.arabidopsis.org) for *Arabidopsis* genes, [UniProt](https://www.uniprot.org) for proteins (Session 3). A GWAS gives **candidates**; only experiments (mutants, swapping alleles between plants) prove that a gene causes the difference. **Which plot, for which question?** | Your question | Plot or number | In R | |---|---|---| | Does one SNP change the trait, and by how much? | boxplot by allele, linear model | `geom_boxplot()`, `lm()` | | Where in the genome are the associations? | Manhattan plot | `geom_point()` + `facet_grid()` | | Can I trust the scan as a whole? | QQ plot, λ | `ppoints()`, `qchisq()` | | Do the plants differ by ancestry? | PCA of the genotypes | `prcomp()` | | Which genes sit under a peak? | the gene table, ±20 kb | `filter()` | ### 7. Worked example: the 1001 Genomes plants - **Flowering time** (1001 Genomes Consortium, 2016): each accession was grown in growth chambers at a constant **10 °C** or **16 °C**, 16 hours of light, **without any winter**, 4 plants per accession. The trait is the number of days from sowing to the first open flower (the mean of the 4 plants). - **Leaf sodium** (Baxter et al., 2010): plants of about 150 of these accessions were grown in the same soil, and the sodium in their leaves measured by mass spectrometry. Too much sodium is toxic for plants: it is the main problem of salty soils. - **Genetic groups**: the 1001 Genomes project grouped the accessions by ancestry (`north_sweden`, `spain`, `relict`, ...); accessions with mixed ancestry are `admixed`. You will use these groups only to look at the results, never in the GWAS itself. ### Practical briefing Blocks A-E in pairs, block F on your own. Block B meets the plants without any statistics: where they come from, and when they flower. Block C tests two SNPs by hand: one that does nothing, and one that looks important but lies. Block D tests all the SNPs, flowering at 16 °C; block E finds the structure with a PCA, corrects the GWAS, and looks up the gene under the best peak. Block F does the same alone for leaf sodium: same code, a different answer. Your submission is the whole rendered workbook: every block counts. ## Your R toolbox New tools today, on top of the Session 4 and Session 5 toolboxes (still valid: see the [Session 4 R toolbox](04_biological_tables_and_visualisation.qmd#your-r-toolbox)). ::: panel-tabset ## A function and apply() A **function** packs a few lines of code under a name, so that you can run them again on other data: ```r snp_p_value <- function(allele) { fit <- lm(flowering ~ allele) summary(fit)$coefficients["allele", "Pr(>|t|)"] } ``` - `allele` is the **argument**: the data you give the function each time you call it; - the last line is what the function returns: here the p-value of the row `allele` in the table of coefficients of `summary(fit)`; - `flowering` is not an argument: the function finds it in your session, so create it first. `apply()` runs a function on every row (or column) of a matrix: ```r apply(genotype_matrix, 1, snp_p_value) # 1 = each row, 2 = each column ``` It returns one value per row: here one p-value per SNP, in the order of the rows. **Common mistake.** `object 'flowering' not found`: run the line that creates `flowering` before `apply()`. ## Matrices: rows, columns, t() | Code | What it does | |---|---| | `genotype_matrix["chr1_920099", ]` | one row (one SNP), all columns (all accessions) | | `genotype_matrix[, 1:5]` | all rows, the first 5 columns | | `dim(genotype_matrix)` | number of rows, then of columns | | `t(genotype_matrix)` | turns the matrix around: rows become columns | | `rownames()`, `colnames()` | the names of the rows and of the columns | `prcomp()` wants the samples (here the accessions) in rows, hence `t()`. ## Manhattan plot ```r ggplot(gwas, aes(x = position / 1e6, y = -log10(p_value))) + geom_point(size = 0.5, alpha = 0.5) + geom_hline(yintercept = -log10(threshold), colour = "red", linetype = "dashed") + facet_grid(~ chromosome, scales = "free_x", space = "free_x") ``` - `position / 1e6`: positions in millions of bases (Mb); - `facet_grid(~ chromosome, ...)`: one panel per chromosome, side by side; `scales = "free_x"` gives each panel its own x axis, `space = "free_x"` makes each panel as wide as its chromosome is long. ## QQ plot and λ ```r qq <- gwas %>% arrange(p_value) %>% mutate(expected = -log10(ppoints(n())), observed = -log10(p_value)) ``` - `ppoints(n())` gives `n` values spread evenly between 0 and 1: the p-values you would expect if no SNP had any effect; - `geom_abline(intercept = 0, slope = 1)` draws the diagonal. The inflation factor: ```r median(qchisq(p_values, df = 1, lower.tail = FALSE)) / qchisq(0.5, df = 1) ``` `qchisq()` turns each p-value back into the test statistic it came from; λ compares the median statistic with the one expected under "no effect". ## A map ```r world <- map_data("world") ggplot() + geom_polygon(data = world, aes(x = long, y = lat, group = group), fill = "grey90", colour = "grey60") + geom_point(data = accessions, aes(x = longitude, y = latitude)) + coord_quickmap(xlim = c(-12, 90), ylim = c(30, 68)) ``` - `map_data("world")` gives the outlines of the countries (from the package maps); - `group = group` draws each country as its own shape; - `coord_quickmap()` keeps the proportions of a map, and `xlim`/`ylim` zoom in. ::: ## Practical (≈1 h 50) ::: callout-tip ## Tested on LEGcompute Everything used today is already installed in RStudio Server's R (tidyverse, maps). Rendering the complete workbook takes about 1.5 to 2 minutes, because of the PCA and the four scans of all the SNPs. ::: ### A. Same habits, new data (15:00-15:15) **A1. The project.** In the **Terminal** tab: ```bash mkdir -p ~/bioinfo_course/session_06/data ~/bioinfo_course/session_06/results cp /home/cmaslard/bioinfo_tools_course/session_06/arabidopsis_accessions.csv ~/bioinfo_course/session_06/data/ cp /home/cmaslard/bioinfo_tools_course/session_06/arabidopsis_genotypes.csv.gz ~/bioinfo_course/session_06/data/ cp /home/cmaslard/bioinfo_tools_course/session_06/arabidopsis_genes.csv ~/bioinfo_course/session_06/data/ cp /home/cmaslard/bioinfo_tools_course/session_06/WORKBOOK.qmd ~/bioinfo_course/session_06/ ls -R ~/bioinfo_course/session_06 ``` The genotype file is compressed (`.gz`). Have a look at its first lines without opening it, still in the Terminal: ```bash zcat ~/bioinfo_course/session_06/data/arabidopsis_genotypes.csv.gz | head -n 4 | cut -c 1-80 ``` Then **File → New Project → Existing Directory**, pick `~/bioinfo_course/session_06`, open `WORKBOOK.qmd`, write your name at the top and click **Render** once. ::: callout-note ## Checkpoint `WORKBOOK.pdf` appears (RStudio opens it), and the first chunk lists the columns `accession_id`, `name`, `country`, `latitude`, `longitude`, `genetic_group`, `flowering_10C`, `flowering_16C`, `leaf_sodium`. The chunks further down show `Error: unexpected input`: those are the `___` you will fill. The `zcat` line prints `snp_id,chromosome,position,acc_88,...` then lines of 0 and 1. ::: `read_csv()` reads the compressed file directly: no need to decompress it. ### B. The accessions, without statistics (15:15-15:35, pairs) **B1. How many, and which?** ```r # how many accessions per genetic group? accessions %>% count(___, sort = TRUE) # how many accessions have a value for each trait? (!is.na() is TRUE when a value is there) accessions %>% summarise(flowering_10C = sum(!is.na(flowering_10C)), flowering_16C = sum(!is.na(flowering_16C)), leaf_sodium = sum(!is.na(___))) ``` ::: callout-note ## Checkpoint Two small tables: the first has one row per genetic group, the second one row with three counts. ::: **B2. Where do they come from?** A map, one point per accession: ```r world <- map_data("world") # the outlines of the countries # one point per accession, coloured by genetic group ggplot() + geom_polygon(data = world, aes(x = long, y = lat, group = group), fill = "grey90", colour = "grey60", linewidth = 0.2) + geom_point(data = accessions, aes(x = longitude, y = latitude, colour = ___), size = 1) + coord_quickmap(xlim = c(-12, 90), ylim = c(30, 68)) + labs(x = "Longitude", y = "Latitude", colour = "Genetic group") ``` ::: callout-note ## Checkpoint A grey map of Europe, North Africa and Asia, with coloured points and a legend of the genetic groups. A grey map with black points means the `colour` is missing in `aes()`. ::: **B3. When do they flower?** The plants grew at a constant 10 °C or 16 °C, without any winter: ```r # flowering at 16 °C, one box per genetic group, the groups sorted by their median ggplot(accessions, aes(x = reorder(genetic_group, flowering_16C, FUN = median, na.rm = TRUE), y = ___)) + geom_boxplot() + coord_flip() + labs(x = NULL, y = "Days to flowering at 16 °C") # the mean flowering time per group, at both temperatures accessions %>% group_by(genetic_group) %>% summarise(mean_16C = mean(flowering_16C, na.rm = TRUE), mean_10C = mean(___)) %>% arrange(desc(mean_16C)) ``` `reorder(genetic_group, flowering_16C, FUN = median, na.rm = TRUE)` sorts the groups by their median flowering time, so that the boxes are in order. **B4. The missing values.** ```r # the accessions without a flowering time at 16 °C: which groups? accessions %>% filter(is.na(___)) %>% count(genetic_group, sort = TRUE) # did they flower early or late at 10 °C? accessions %>% group_by(missing_16C = is.na(flowering_16C)) %>% summarise(n = n(), mean_10C = mean(flowering_10C)) ``` ### 15:35-15:45: Short break ### C. One SNP at a time (15:45-16:10, pairs) **C1. The genotypes.** ```r # one row per SNP (as in a VCF file): its id, chromosome, position, then one column per accession genotypes <- read_csv("data/arabidopsis_genotypes.csv.gz") dim(genotypes) genotypes[1:5, 1:8] # the first 5 SNPs, the first 5 accessions # the 0/1 only, as a matrix: one row per SNP, one column per accession genotype_matrix <- genotypes %>% select(starts_with(___)) %>% as.matrix() rownames(genotype_matrix) <- genotypes$___ dim(genotype_matrix) # are the accessions in the same order as the rows of the table `accessions`? all(colnames(genotype_matrix) == paste0("acc_", accessions$accession_id)) ``` ::: callout-note ## Checkpoint `dim()` prints two numbers twice: rows, then columns (the table has 3 more columns than the matrix). `genotypes[1:5, 1:8]` shows `snp_id`, `chromosome`, `position` and 0/1 values. The last line must print `TRUE`: if it prints `FALSE`, stop and call the instructor. ::: **C2. A first SNP.** `genotype_matrix["chr1_920099", ]` takes one row of the matrix: the allele of every accession at this SNP, in the order of `accessions`. ```r # a first SNP: its allele in each accession, next to the flowering time one_snp <- accessions %>% select(accession_id, genetic_group, flowering_16C) %>% mutate(allele = genotype_matrix["chr1_920099", ]) one_snp %>% count(___) ggplot(one_snp, aes(x = factor(allele), y = flowering_16C)) + geom_boxplot(outlier.shape = NA) + geom_jitter(width = 0.2, alpha = 0.3) + labs(x = "Allele at chr1_920099", y = "Days to flowering at 16 °C") # the linear model of Session 4: does the allele change the flowering time? summary(lm(flowering_16C ~ ___, data = one_snp)) ``` **C3. A second SNP**, the same way, then coloured by genetic group: ```r # a second SNP, the same way trap_snp <- accessions %>% select(accession_id, genetic_group, flowering_16C) %>% mutate(allele = genotype_matrix[___, ]) summary(lm(flowering_16C ~ allele, data = trap_snp)) # the same plot, the points coloured by genetic group ggplot(trap_snp, aes(x = factor(allele), y = flowering_16C)) + geom_boxplot(outlier.shape = NA) + geom_jitter(aes(colour = ___), width = 0.2, size = 1) + labs(x = "Allele at chr2_13280384", y = "Days to flowering at 16 °C") # the share of accessions with allele 1, per genetic group (the mean of 0s and 1s) trap_snp %>% group_by(genetic_group) %>% summarise(share_allele_1 = mean(___)) %>% arrange(desc(share_allele_1)) ``` `mean()` of 0s and 1s is the share of 1s: `mean(c(0, 1, 1, 1))` is 0.75. **C4. The same model, with the genetic group.** ```r # the same model, with the genetic group added (as heat next to water in Session 4) summary(lm(flowering_16C ~ allele + ___, data = trap_snp)) ``` ::: callout-note ## Checkpoint C2 to C4 each print a `lm()` summary with a row `allele`. C4 has many more rows: one per genetic group (compared with `admixed`, the first in alphabetical order). ::: ### D. Guided: one linear model per SNP (16:10-16:35, pairs) **D1. All the SNPs.** The model of C2, packed into a function, then run on every row of the matrix: ```r # the trait, in the same order as the columns of genotype_matrix flowering <- accessions$flowering_16C # the p-value of one SNP: the lm() of C2, written as a small function snp_p_value <- function(allele) { fit <- lm(flowering ~ allele) summary(fit)$coefficients["allele", "Pr(>|t|)"] } # check it on the SNP of C3: the same p-value as there? snp_p_value(genotype_matrix[___, ]) # now every SNP: apply() runs the function on each row (1) of the matrix (about 20 s) p_naive <- apply(genotype_matrix, ___, snp_p_value) # one row per SNP: its position and its p-value gwas_naive <- genotypes %>% select(snp_id, chromosome, position) %>% mutate(p_value = ___) gwas_naive %>% arrange(p_value) # how many p-values below 0.05? sum(gwas_naive$p_value < ___) ``` ::: callout-note ## Checkpoint The first line prints one p-value; `apply()` then works for about 20 seconds; the table `gwas_naive` has 4 columns (`snp_id`, `chromosome`, `position`, `p_value`) and one row per SNP. `object 'flowering' not found` means that the first line was not run. ::: **D2. The threshold and the Manhattan plot.** ```r # Bonferroni: 5% divided by the number of tests threshold <- 0.05 / ___ threshold -log10(threshold) sum(gwas_naive$p_value < ___) # the Manhattan plot: one point per SNP, one panel per chromosome ggplot(gwas_naive, aes(x = position / 1e6, y = -log10(___))) + geom_point(size = 0.5, alpha = 0.5) + geom_hline(yintercept = -log10(threshold), colour = "red", linetype = "dashed") + facet_grid(~ chromosome, scales = "free_x", space = "free_x") + labs(title = "Flowering at 16 °C, naive GWAS", x = "Position (Mb)", y = "-log10(p)") ``` **D3. The QQ plot and λ.** ```r # expected: the p-values if no SNP had any effect (spread evenly between 0 and 1) qq_naive <- gwas_naive %>% arrange(p_value) %>% mutate(expected = -log10(ppoints(n())), observed = -log10(___)) ggplot(qq_naive, aes(x = ___, y = ___)) + geom_point(size = 0.5) + geom_abline(intercept = 0, slope = 1, colour = "red") + labs(title = "QQ plot, naive GWAS", x = "Expected -log10(p)", y = "Observed -log10(p)") # the inflation factor lambda: 1 means no inflation lambda_naive <- median(qchisq(gwas_naive$p_value, df = 1, lower.tail = FALSE)) / qchisq(0.5, df = 1) lambda_naive ``` ::: callout-note ## Checkpoint D2: one very small number, the same on the −log10 scale, a count, and a plot with five panels (the five chromosomes) and a red dashed line. D3: a plot with a red diagonal, and one number. ::: ### E. Guided: the population structure, then the corrected GWAS (16:35-17:10, pairs) **E1. The PCA of the genotypes.** The same `prcomp()` as in Session 5, on all the accessions and all the SNPs (about 30 seconds): ```r # the PCA of Session 5: accessions in rows (t() turns the matrix), SNPs in columns (about 30 s) genotype_pca <- prcomp(t(___), scale. = ___) pca_variance <- round(100 * summary(genotype_pca)$importance[2, 1:10], 1) pca_variance pca_scores <- as_tibble(genotype_pca$x[, 1:2]) %>% bind_cols(select(accessions, genetic_group, latitude, longitude)) ggplot(pca_scores, aes(x = PC1, y = PC2, colour = ___)) + geom_point(size = 1) + labs(x = paste0("PC1 (", pca_variance[1], "%)"), y = paste0("PC2 (", pca_variance[___], "%)"), colour = "Genetic group") # does PC1 follow geography? cor(pca_scores$PC1, pca_scores$___, use = "complete.obs") cor(pca_scores$PC1, pca_scores$longitude, use = "complete.obs") ``` ::: callout-note ## Checkpoint Ten percentages, a plot with one point per accession coloured by genetic group, and two correlations. On another computer the plot may be flipped left-right or up-down: the direction of an axis is arbitrary. ::: **E2. How much of flowering follows ancestry?** ```r # the first 10 PCs, as a matrix: one row per accession, one column per PC pcs <- genotype_pca$x[, 1:___] # how much of the flowering time do they explain? (R², the share of variation of Session 4) summary(lm(flowering ~ ___))$r.squared ``` **E3. The corrected GWAS.** The function of D1, with the 10 PCs added to the model: ```r # the same function, with the 10 PCs in the model (as genetic_group in C4) snp_p_value_pc <- function(allele) { fit <- lm(flowering ~ allele + ___) summary(fit)$coefficients["allele", "Pr(>|t|)"] } p_pc <- apply(genotype_matrix, 1, ___) # about 30 s gwas_pc <- genotypes %>% select(snp_id, chromosome, position) %>% mutate(p_value = ___) sum(gwas_pc$p_value < ___) lambda_pc <- median(qchisq(gwas_pc$p_value, df = 1, lower.tail = FALSE)) / qchisq(0.5, df = 1) lambda_pc # the Manhattan plot (as in D2) ggplot(gwas_pc, aes(x = position / 1e6, y = -log10(___))) + geom_point(size = 0.5, alpha = 0.5) + geom_hline(yintercept = -log10(___), colour = "red", linetype = "dashed") + facet_grid(~ ___, scales = "free_x", space = "free_x") + labs(title = "Flowering at 16 °C, GWAS with 10 PCs", x = "Position (Mb)", y = "-log10(p)") # the QQ plots of both scans, on one plot qq_both <- bind_rows(naive = gwas_naive, corrected = gwas_pc, .id = "model") %>% group_by(model) %>% arrange(p_value, .by_group = TRUE) %>% mutate(expected = -log10(ppoints(n())), observed = -log10(p_value)) ggplot(qq_both, aes(x = expected, y = observed, colour = ___)) + geom_point(size = 0.5) + geom_abline(intercept = 0, slope = 1) + labs(title = "QQ plot, before and after the correction", x = "Expected -log10(p)", y = "Observed -log10(p)") ``` **E4. The best peak, and its genes.** ```r # the SNP of C3, before and after the correction gwas_naive %>% filter(snp_id == "chr2_13280384") gwas_pc %>% filter(snp_id == ___) # the 10 best SNPs after the correction gwas_pc %>% arrange(p_value) %>% head(10) # the genes within 20 kb (20,000 bp) of the best SNP genes <- read_csv("data/arabidopsis_genes.csv") top_snp <- gwas_pc %>% slice_min(___, n = 1) top_snp genes %>% filter(chromosome == top_snp$___, end > top_snp$position - ___, start < top_snp$position + 20000) # the same SNP in the naive scan: its p-value, and how many SNPs had a smaller one naive_p_top <- gwas_naive %>% filter(snp_id == top_snp$___) %>% pull(p_value) naive_p_top sum(gwas_naive$p_value < naive_p_top) ``` **E4, continued: the gene names on the plot.** Nothing to fill: after your answer to E4, run this chunk. It draws the Manhattan plot of E3 again, with the best SNP in blue and the short names of the genes within 20 kb of it. ```r # the short names of the genes within 20 kb of the best SNP (the same genes as above) top_genes <- genes %>% filter(chromosome == top_snp$chromosome, end > top_snp$position - 20000, start < top_snp$position + 20000, !is.na(symbol)) %>% mutate(symbol = str_replace_all(symbol, "%3B", ";")) # the file writes ";" as %3B top_label <- top_snp %>% mutate(genes = paste(top_genes$symbol, collapse = "\n")) # one name per line # the Manhattan plot of E3, with the best SNP in blue and the names of its genes ggplot(gwas_pc, aes(x = position / 1e6, y = -log10(p_value))) + geom_point(size = 0.5, alpha = 0.5) + geom_hline(yintercept = -log10(threshold), colour = "red", linetype = "dashed") + geom_point(data = top_label, colour = "blue", size = 2) + geom_label(data = top_label, aes(label = genes), hjust = 0, vjust = 1, nudge_x = 1, size = 2.5) + facet_grid(~ chromosome, scales = "free_x", space = "free_x") + labs(title = "The genes near the best SNP (10 PCs)", x = "Position (Mb)", y = "-log10(p)") ``` ::: callout-note ## Checkpoint E3: a count, one number (λ), a Manhattan plot titled "GWAS with 10 PCs", and a QQ plot with two colours. E4: tables of 4 columns, then a table of genes with `gene_id`, `symbol`, `chromosome`, `start`, `end`, `description`. A `symbol` shown as `NA` means the gene has no short name. E4, continued: the same Manhattan plot, with one blue point and a box of gene names next to it. ::: Then **E5**: summarise the flowering GWAS in a few sentences, with your own numbers. ### F. On your own: leaf sodium The same reasoning, on the **leaf sodium** of about 150 accessions. The workbook gives you the skeleton, with more `___` than in D and E: by now, most of them are the choices you made there, so reuse your own code. - **F1**: the distribution of sodium, how much the 3 first PCs explain, and the naive scan with its λ. `lm()` leaves out the accessions without a sodium value (`NA`) by itself. With 147 plants, 3 PCs are enough: check λ. - **F2**: the corrected scan and its Manhattan plot, saved with `ggsave()` and inserted into your report, as in Session 5, with your own caption and the label `#fig-sodium`. - **F3**: the genes under the best SNP, and the sodium of the plants by allele. - **F4**, comparing the two GWAS, and **F5**, one limitation. ### Wrap-up (17:10-17:15) Click **Render** (about 1.5 to 2 minutes for the whole workbook). Check that your PDF shows all the plots of blocks B to E. Block F follows the same logic as D and E. Submit `WORKBOOK.pdf` on Moodle **by Sunday evening (1 November)**. If you see errors in your PDF, or if it does not render, submit in addition your `WORKBOOK.qmd` and your `results/` folder (your saved plot). ::: callout-note ## How this is graded (10 points) | Criterion | Points | What it looks for | |---|---:|---| | Reproducible workflow | 2 | The workbook renders (or the code is complete and in order), inside a clean project. | | Working analysis | 3 | All the code of blocks A to F ran, and all your answers use the right numbers. | | Required result | 2 | The corrected Manhattan plot of leaf sodium (F2): saved, inserted with a caption and `#fig-sodium`, cited in F4. | | Biological interpretation | 2 | Your answers on the biology, block by block, and especially E5 and F4, with your own numbers. | | Reflection | 1 | Your answers on the methods (the order of the accessions, the threshold, the QQ plot, λ, the correction) and your limitation (F5). | Every question counts, except the exit questions. ::: ## Checkpoint You are done when all of the following are true: - you can say what a SNP is and what 0 and 1 mean in the genotype table; - you tested one SNP with `lm()`, then all of them with a function and `apply()`; - you can explain the Bonferroni threshold, and read a Manhattan plot; - you can read a QQ plot and λ, and say when a scan cannot be trusted; - you can explain, with the SNP of C3, how population structure creates false associations; - you ran a PCA of the genotypes, and used its components to correct the GWAS; - you found the genes under a peak, and can say why a peak is not a proof; - you can explain why the structure misleads the GWAS of flowering time much more than that of leaf sodium. ## Optional challenge For your own interest, not graded. **A. How many PCs?** Redo E3 with 3, then 20 PCs instead of 10. How do λ and the number of SNPs below the threshold change? Is the best SNP still the same? **B. *FRIGIDA*.** *FRI* (chromosome 4, around 0.27 Mb) is one of the best-known flowering genes of *Arabidopsis*. Find the best SNP within 20 kb of it in the naive scan and in the corrected scan. What happened to it, and why could a correction for ancestry also hide a real gene? **C. Flowering at 10 °C.** Redo D and E with `flowering_10C`. Is the best peak in the same place as at 16 °C? **D. A prettier Manhattan plot.** Colour the points by chromosome, and label the best SNP of each tower with `geom_text_repel()` (package ggrepel, Session 5). ## Going further - Korte A. and Farlow A. (2013). [The advantages and limitations of trait analysis with GWAS: a review](https://doi.org/10.1186/1746-4811-9-29). *Plant Methods*. The best short review for plant biologists. - 1001 Genomes Consortium (2016). [1,135 genomes reveal the global pattern of polymorphism in *Arabidopsis thaliana*](https://doi.org/10.1016/j.cell.2016.05.063). *Cell*. Where today's SNPs and flowering times come from. - [AraGWAS Catalog](https://aragwas.1001genomes.org): published *Arabidopsis* GWAS, including flowering time, run with mixed models: compare with your peaks. - [AraPheno](https://arapheno.1001genomes.org): hundreds of other traits measured on the same accessions. ## Before the next session Session 7 moves from testing to **predicting**: deep learning. You will meet again the idea behind the cross-validation of Session 5: a model is only judged on data it has never seen.