4 Genome-Wide Association: a Thrombophilia Case

In this lab we draft a genome-wide association (GWA) study. Genome-wide association is the overall name for a technique that seeks positions in the genome associated with a given phenotype. We carry out this laboratory with the maintained gaston package (Hervé Perdry & Claire Dandine-Roulland), which fits the linear mixed models needed to account for family relatedness.¹

By the end of this laboratory you will be able to:

load a real family-based GWA dataset in PLINK/bed format and explore it
look up SNP markers in dbSNP
reason about familial structure and kinship, and visualise it
estimate trait heritability and run a mixed-model genome-wide scan that accounts for relatedness
read a Manhattan and QQ plot, and identify and interpret the top associated SNP

Lab at a glance

Goal: run a family-based GWAS with a linear mixed model and reason about relatedness and confidentiality. Prerequisites: Genetic Association. Time: ~2 h. Packages: gaston. Data: the GAIT dataset in gaston/bed format (see below). Deliverable: a lab report.

For this session we use real data from the GAIT project (Genetic Analysis of Idiopathic Thrombophilia), kindly provided by Dr. Soria (Unit of Genomics and Complex Diseases, Institut de Recerca de l’Hospital de la Santa Creu i Sant Pau). GAIT is an NIH-funded project studying strong linkage signals for QTLs influencing quantitative risk factors for thrombosis, detected in the 21 extended families that comprise the GAIT sample.

Confidential data

This dataset is not public. The information it contains is privileged, confidential, and protected from disclosure. Unauthorised use, reproduction, or distribution of the dataset — or any portion of it, or any result derived from it — is strictly forbidden, may result in severe civil and criminal penalties, and will be prosecuted to the maximum extent possible under the law. If you received this data in error, notify us immediately (alexandre.perera@upc.edu).

The genotypes are stored in the widely used PLINK bed format, which gaston reads by slices. A bed dataset is three files: .bed (genotypes), .fam (individuals), and .bim (marker map). If you only have the legacy LabGWAS1.Rdata, convert it once with the helper script (Rscript scripts/make-gaston-bed.R), which writes objQ2.bed/.bim/.fam and objQ2.phe into data/.

Code

library(gaston)
geno  <- read.bed.matrix("data/objQ2")

Reading data/objQ2.rds 
Reading data/objQ2.bed

Code

pheno <- read.table("data/objQ2.phe", header = TRUE)

geno is a bed.matrix: two annotation tables, @ped (individuals) and @snps (markers), plus a pointer to the packed genotypes. pheno holds the phenotype table.

4.1 Data exploration

Code

geno

A bed.matrix with 377 individuals and 266966 markers.
snps stats are set
ped stats are set

The @ped table counts, per individual, the genotypes equal to 0/1/2 and missing (N0, N1, N2, NAs), the call rate and heterozygosity; @snps holds the alternate allele frequency p, its expected genotype mu = 2p, and the genotype standard deviation sigma. Under Hardy–Weinberg equilibrium sigma should track $\sqrt{2p(1-p)}$ (Figure 4.1):

Code

plot(geno@snps$p, geno@snps$sigma, xlim = c(0, 1), pch = ".",
     xlab = "allele frequency p", ylab = expression(sigma))
t <- seq(0, 1, length = 101)
lines(t, sqrt(2 * t * (1 - t)), col = "#c0432f", lwd = 2)

Scatter of sigma versus allele frequency p, hugging the sqrt(2p(1-p)) Hardy-Weinberg curve. — Figure 4.1: Genotype standard deviation against alternate-allele frequency; the red curve is the Hardy–Weinberg expectation.

✅ Exercise 4.1 — Data exploration

What is the number of individuals in the study? And the number of families (pheno$FAM)?
How many members were recruited from each family?
Does the $p$–$\sigma$ cloud follow the H–W expectation? Where would you expect departures?
Build a table of the Smoking (SMOKIN) vs Gender (sex) data.

Solution 4.1

377 individuals across 22 families (mean ~17 per family).
The $p$–$\sigma$ cloud broadly hugs the $\sqrt{2p(1-p)}$ curve; departures appear at the highest allele frequencies (and from genotyping error or non-random mating in these large families).
Smoking by gender: table(pheno$SMOKIN, pheno$sex) gives (row = SMOKIN 1/2, col = sex 0/1) 142/89 and 62/82.

4.1.1 Genotype data

The genotypes cast to a numeric 0/1/2 matrix with as.matrix(); the rs codes are the column names. The first few markers:

Code

as.matrix(geno)[1:5, 1:5]

     rs3934834 rs3737728 rs6687776 rs9651273 rs4970405
1203         0         1         0         1         0
1204         0         0         0         1         0
1208         0         1         0         1         0
1209         0         1         0         1         0
1210         0         1         0         2         0

Look up a marker on the dbSNP website (https://www.ncbi.nlm.nih.gov/snp/) — a public catalogue of Single Nucleotide Polymorphisms. Marker names are rs accessions (Reference SNP records).

✅ Exercise 4.2 — Markers

On which chromosome and position is the first SNP located (geno@snps[1, ])?

Solution 4.2

geno@snps[1, ] is rs3934834 on chromosome 1, position 995669 (a known 1p36 SNP you can confirm in dbSNP).

4.2 Familial structure

Each individual carries family (FAM) and parent (FA, MO) information. Most individuals are related, and this kinship must be accounted for when computing association statistics — otherwise relatedness inflates the test (recall the genomic-inflation $\lambda$ from the Genetic Association lab).

We summarise relatedness with the Genetic Relationship Matrix (GRM), the standardised whole-genome analogue of the Identity-By-State idea: for individuals $i, j$ over $N$ SNPs with allele frequency $p_k$,

\[ K_{i,j} = \frac{1}{N} \sum_{k=1}^{N} \frac{(x_{i,k} - 2p_k)(x_{j,k} - 2p_k)}{2 p_k (1 - p_k)} \tag{4.1}\]

Code

K <- GRM(geno)
K[1:5, 1:5]

             1203         1204         1208          1209          1210
1203  0.988955021 -0.005105684  0.486438155 -0.0178462826  0.3734436631
1204 -0.005105684  1.000756979 -0.008204872  0.0072001172 -0.0163245704
1208  0.486438155 -0.008204872  0.990921259 -0.0019851427  0.4473317862
1209 -0.017846283  0.007200117 -0.001985143  1.0205782652 -0.0006051082
1210  0.373443663 -0.016324570  0.447331786 -0.0006051082  0.9902690649

A hierarchical clustering of $1 - K/\max K$ recovers the family blocks (Figure 4.2):

Code

d  <- as.dist(1 - K / max(K))
hc <- hclust(d, method = "single")
plot(hc, labels = geno@ped$famid, cex = 0.4, main = "", xlab = "", sub = "")

Dendrogram of individuals clustering into family groups, leaves labelled by family id. — Figure 4.2: Hierarchical clustering from the genetic relationship matrix, labelled by family — members of a pedigree group together.

✅ Exercise 4.3 — Kinship

How many unrelated individuals (founders) are in the dataset (FA/MO missing)?
Why is the GRM computed from autosomal SNPs only?
Comment on Figure 4.2: do the clusters match the families?

Solution 4.3

94 founders (no parents in the data).
Autosomes are diploid in everyone and largely independent of sex, so their kinship model is valid for both sexes; sex chromosomes follow different (hemizygous) inheritance and would bias relatedness estimates.
The dendrogram blocks group individuals of the same famid together — the GRM recovers the families.

4.3 Association

The target phenotype XXX is a critical, highly heritable clinical trait related to hemostasis. Note that XXX is a deliberately anonymised label: because the dataset is under confidentiality the real trait name is replaced by a placeholder, so treat it as a generic quantitative phenotype and focus on the method (kinship accounting), not the specific analyte. We first fit the polygenic model — a linear mixed model with a random effect structured by the GRM — to estimate how much of the trait variance is genetic (heritability, $h^2$). We restrict to individuals with complete phenotype and covariates:

Code

Y   <- pheno$XXX
Xc  <- cbind(AGE = pheno$AGE, sex = pheno$sex)
ok  <- complete.cases(Y, Xc)
fit <- lmm.aireml(Y[ok], X = cbind(1, Xc[ok, ]), K = K[ok, ok], verbose = FALSE)
h2  <- fit$tau / (fit$tau + fit$sigma2)
round(c(genetic_var = fit$tau, resid_var = fit$sigma2, heritability = h2), 3)

 genetic_var    resid_var heritability 
     900.390      196.472        0.821

The heritability is high — most of the variance in XXX is genetic, which is why a relatedness-aware model matters. Now the genome-wide mixed-model score test (gaston’s equivalent of GenABEL’s mmscore): each SNP is tested for association with XXX, adjusting for age, sex, and the full kinship:

Code

assoc <- association.test(geno[ok, ], Y = Y[ok], X = Xc[ok, ],
                          method = "lmm", K = K[ok, ok], test = "score")

[Iteration 1] theta = 533.503 537.379
[Iteration 1] log L = -1434.21
[Iteration 1] AI-REML update
[Iteration 1] ||gradient|| = 0.0351275
[Iteration 2] theta = 133.062 887.195
[Iteration 2] log L = -1428.5
[Iteration 2] AI-REML update
[Iteration 2] ||gradient|| = 0.0477282
[Iteration 3] theta =  182.05 911.125
[Iteration 3] log L = -1427.16
[Iteration 3] AI-REML update
[Iteration 3] ||gradient|| = 0.00529741
[Iteration 4] theta = 195.403 901.948
[Iteration 4] log L = -1427.13
[Iteration 4] AI-REML update
[Iteration 4] ||gradient|| = 0.000214138
[Iteration 5] theta = 196.418 900.487
[Iteration 5] log L = -1427.13
[Iteration 5] AI-REML update
[Iteration 5] ||gradient|| = 8.26671e-06

Code

assoc <- assoc[!is.na(assoc$p), ]

A Manhattan plot shows $-\log_{10} p$ along the genome (Figure 4.3):

Code

manhattan(assoc, main = "XXX association (mixed model)")
abline(h = -log10(5e-8), col = "#c0432f", lty = 2)

Manhattan plot with one clearly significant peak rising far above the genome-wide significance threshold on chromosome 5. — Figure 4.3: Manhattan plot of the mixed-model association to the XXX phenotype. The dashed line is genome-wide significance (5e-8).

The scan is well-calibrated because relatedness is modelled (Figure 4.4):

Code

qqplot.pvalues(assoc$p)

QQ plot with points on the diagonal except a short tail of genuine associations. — Figure 4.4: QQ plot of the mixed-model association p-values.

The top associated SNPs:

Code

head(assoc[order(assoc$p), c("chr", "pos", "id", "A1", "A2", "freqA2", "p")], 10)

✅ Exercise 4.4 — Top hit

Locate the SNP with the smallest p-value — its rs id, chromosome, and position. Look it up in dbSNP.
Which gene lies at that locus, and what is its role in coagulation/thrombosis?
Compare the mixed-model result with a naïve test that ignores kinship (association.test(geno[ok,], Y[ok], X = Xc[ok,], method = "lm")). What changes, and why does ignoring relatedness mislead?

Solution 4.4

The top hit is rs2731672 (chr 5, ~176.8 Mb), $p \approx 2.2\times10^{-27}$, an intergenic/regulatory SNP near F12 (coagulation factor XII, the contact-activation protease of the intrinsic pathway — directly relevant to hemostasis/thrombosis).
The modelled heritability is high ($h^2 \approx 0.82$). A naïve per-SNP linear test ignores the family correlation, so observations are not independent → p-values are anti-conservative and inflated. The mixed model uses the GRM to absorb the polygenic/family background, leaving a well-calibrated scan (the QQ plot hugs the diagonal except the genuine tail).

Ethics — genomic data privacy

Genotype data are inherently identifying: a few dozen SNPs can re-identify an individual, and a genome reveals information about relatives who never consented. That is why the GAIT set is access-controlled and why real GWAS data live in managed-access repositories (dbGaP, EGA) under data-use agreements. Share summary statistics — never individual genotypes — and follow the study’s consent and governance terms.

4.4 Further reading

Zhou & Stephens (2012), Genetics — the polygenic linear mixed model (lmm.aireml) implemented here.
gaston vignette “linear mixed model” — association.test(..., method = "lmm").
Soria et al., the GAIT project papers, for the thrombophilia phenotype background.

Worked solutions for the exercises are in the collapsed callouts below each exercise.

4.5 Report

Write a laboratory report following the structure of this text. It should contain an introduction, the goals of each section, and a description of the results of each part. Include any graphics relevant to making your work understandable. The lab should be delivered within two weeks.

Earlier editions of this lab used GenABEL (polygenic + mmscore); GenABEL has since been archived from CRAN. gaston reproduces the same family-aware analysis — lmm.aireml estimates the polygenic variance components and association.test(..., method = "lmm") is the mixed-model score test — and is what we use here.↩︎

--- title: "Genome-Wide Association: a Thrombophilia Case" --- In this lab we draft a genome-wide association (GWA) study. Genome-wide association is the overall name for a technique that seeks positions in the genome *associated* with a given phenotype. We carry out this laboratory with the maintained **gaston** package (Hervé Perdry & Claire Dandine-Roulland), which fits the linear mixed models needed to account for family relatedness.[^gaston-note] [^gaston-note]: Earlier editions of this lab used **GenABEL** (`polygenic` + `mmscore`); GenABEL has since been archived from CRAN. `gaston` reproduces the same family-aware analysis — `lmm.aireml` estimates the polygenic variance components and `association.test(..., method = "lmm")` is the mixed-model score test — and is what we use here. By the end of this laboratory you will be able to: - load a real family-based GWA dataset in PLINK/bed format and explore it - look up SNP markers in dbSNP - reason about familial structure and kinship, and visualise it - estimate trait heritability and run a **mixed-model** genome-wide scan that accounts for relatedness - read a Manhattan and QQ plot, and identify and interpret the top associated SNP ```{r} #| label: gait-setup #| include: false source("R/setup.R") ``` ::: {.callout-note appearance="minimal" title="Lab at a glance"} **Goal:** run a family-based GWAS with a linear mixed model and reason about relatedness and confidentiality. **Prerequisites:** [Genetic Association](association.qmd). **Time:** ~2 h. **Packages:** `gaston`. **Data:** the GAIT dataset in gaston/bed format (see below). **Deliverable:** a lab report. ::: For this session we use real data from the **GAIT** project (Genetic Analysis of Idiopathic Thrombophilia), kindly provided by Dr. Soria (Unit of Genomics and Complex Diseases, Institut de Recerca de l'Hospital de la Santa Creu i Sant Pau). GAIT is an NIH-funded project studying strong linkage signals for QTLs influencing quantitative risk factors for thrombosis, detected in the 21 extended families that comprise the GAIT sample. ::: {.callout-important title="Confidential data"} This dataset is **not public**. The information it contains is privileged, confidential, and protected from disclosure. Unauthorised use, reproduction, or distribution of the dataset — or any portion of it, or any result derived from it — is strictly forbidden, may result in severe civil and criminal penalties, and will be prosecuted to the maximum extent possible under the law. If you received this data in error, notify us immediately (<alexandre.perera@upc.edu>). ::: The genotypes are stored in the widely used **PLINK bed** format, which `gaston` reads by slices. A bed dataset is three files: `.bed` (genotypes), `.fam` (individuals), and `.bim` (marker map). If you only have the legacy `LabGWAS1.Rdata`, convert it once with the helper script (`Rscript scripts/make-gaston-bed.R`), which writes `objQ2.bed/.bim/.fam` and `objQ2.phe` into `data/`. ```{r} #| label: gait-load #| message: false library(gaston) geno <- read.bed.matrix("data/objQ2") pheno <- read.table("data/objQ2.phe", header = TRUE) ``` `geno` is a **bed.matrix**: two annotation tables, `@ped` (individuals) and `@snps` (markers), plus a pointer to the packed genotypes. `pheno` holds the phenotype table. ## Data exploration {#sec-gait-data-exploration} ```{r} #| label: gait-show geno ``` The `@ped` table counts, per individual, the genotypes equal to 0/1/2 and missing (`N0, N1, N2, NAs`), the call rate and heterozygosity; `@snps` holds the alternate allele frequency `p`, its expected genotype `mu = 2p`, and the genotype standard deviation `sigma`. Under Hardy–Weinberg equilibrium `sigma` should track $\sqrt{2p(1-p)}$ (@fig-gait-hwe): ```{r} #| label: fig-gait-hwe #| fig-cap: "Genotype standard deviation against alternate-allele frequency; the red curve is the Hardy–Weinberg expectation." #| fig-alt: "Scatter of sigma versus allele frequency p, hugging the sqrt(2p(1-p)) Hardy-Weinberg curve." #| fig-width: 5.5 #| fig-height: 4 plot(geno@snps$p, geno@snps$sigma, xlim = c(0, 1), pch = ".", xlab = "allele frequency p", ylab = expression(sigma)) t <- seq(0, 1, length = 101) lines(t, sqrt(2 * t * (1 - t)), col = "#c0432f", lwd = 2) ``` ::: {.callout-tip title="✅ Exercise 4.1 — Data exploration"} 1. What is the number of individuals in the study? And the number of families (`pheno$FAM`)? 2. How many members were recruited from each family? 3. Does the $p$–$\sigma$ cloud follow the H–W expectation? Where would you expect departures? 4. Build a table of the Smoking (`SMOKIN`) vs Gender (`sex`) data. ::: ::: {.content-hidden when-format="pdf" unless-profile="solutions"} ::: {.callout-tip collapse="true" title="Solution 4.1"} - **377 individuals** across **22 families** (mean ~17 per family). - The $p$–$\sigma$ cloud broadly hugs the $\sqrt{2p(1-p)}$ curve; departures appear at the highest allele frequencies (and from genotyping error or non-random mating in these large families). - Smoking by gender: `table(pheno$SMOKIN, pheno$sex)` gives (row = SMOKIN 1/2, col = sex 0/1) `142/89` and `62/82`. ::: ::: ### Genotype data {#sec-gait-genotype-data} The genotypes cast to a numeric 0/1/2 matrix with `as.matrix()`; the `rs` codes are the column names. The first few markers: ```{r} #| label: gait-first-markers as.matrix(geno)[1:5, 1:5] ``` Look up a marker on the dbSNP website (<https://www.ncbi.nlm.nih.gov/snp/>) — a public catalogue of Single Nucleotide Polymorphisms. Marker names are `rs` accessions (Reference SNP records). ::: {.callout-tip title="✅ Exercise 4.2 — Markers"} On which chromosome and position is the first SNP located (`geno@snps[1, ]`)? ::: ::: {.content-hidden when-format="pdf" unless-profile="solutions"} ::: {.callout-tip collapse="true" title="Solution 4.2"} `geno@snps[1, ]` is **rs3934834** on **chromosome 1**, position **995669** (a known 1p36 SNP you can confirm in dbSNP). ::: ::: ## Familial structure {#sec-familial-structure} Each individual carries family (`FAM`) and parent (`FA`, `MO`) information. Most individuals are related, and this kinship must be accounted for when computing association statistics — otherwise relatedness inflates the test (recall the genomic-inflation $\lambda$ from the [Genetic Association](association.qmd#sec-qq-initial) lab). We summarise relatedness with the **Genetic Relationship Matrix** (GRM), the standardised whole-genome analogue of the Identity-By-State idea: for individuals $i, j$ over $N$ SNPs with allele frequency $p_k$, $$ K_{i,j} = \frac{1}{N} \sum_{k=1}^{N} \frac{(x_{i,k} - 2p_k)(x_{j,k} - 2p_k)}{2 p_k (1 - p_k)} $$ {#eq-grm} ```{r} #| label: gait-grm K <- GRM(geno) K[1:5, 1:5] ``` A hierarchical clustering of $1 - K/\max K$ recovers the family blocks (@fig-gait-ibs): ```{r} #| label: fig-gait-ibs #| fig-cap: "Hierarchical clustering from the genetic relationship matrix, labelled by family — members of a pedigree group together." #| fig-alt: "Dendrogram of individuals clustering into family groups, leaves labelled by family id." #| fig-width: 9 #| fig-height: 4 d <- as.dist(1 - K / max(K)) hc <- hclust(d, method = "single") plot(hc, labels = geno@ped$famid, cex = 0.4, main = "", xlab = "", sub = "") ``` ::: {.callout-tip title="✅ Exercise 4.3 — Kinship"} 1. How many unrelated individuals (founders) are in the dataset (`FA`/`MO` missing)? 2. Why is the GRM computed from autosomal SNPs only? 3. Comment on @fig-gait-ibs: do the clusters match the families? ::: ::: {.content-hidden when-format="pdf" unless-profile="solutions"} ::: {.callout-tip collapse="true" title="Solution 4.3"} - **94 founders** (no parents in the data). - Autosomes are diploid in everyone and largely independent of sex, so their kinship model is valid for both sexes; sex chromosomes follow different (hemizygous) inheritance and would bias relatedness estimates. - The dendrogram blocks group individuals of the same `famid` together — the GRM recovers the families. ::: ::: ## Association {#sec-gait-association} The target phenotype `XXX` is a critical, highly heritable clinical trait related to hemostasis. Note that `XXX` is a deliberately **anonymised label**: because the dataset is under confidentiality the real trait name is replaced by a placeholder, so treat it as a generic quantitative phenotype and focus on the *method* (kinship accounting), not the specific analyte. We first fit the **polygenic model** — a linear mixed model with a random effect structured by the GRM — to estimate how much of the trait variance is genetic (heritability, $h^2$). We restrict to individuals with complete phenotype and covariates: ```{r} #| label: gait-polygenic Y <- pheno$XXX Xc <- cbind(AGE = pheno$AGE, sex = pheno$sex) ok <- complete.cases(Y, Xc) fit <- lmm.aireml(Y[ok], X = cbind(1, Xc[ok, ]), K = K[ok, ok], verbose = FALSE) h2 <- fit$tau / (fit$tau + fit$sigma2) round(c(genetic_var = fit$tau, resid_var = fit$sigma2, heritability = h2), 3) ``` The heritability is high — most of the variance in `XXX` is genetic, which is why a relatedness-aware model matters. Now the genome-wide **mixed-model score test** (`gaston`'s equivalent of GenABEL's `mmscore`): each SNP is tested for association with `XXX`, adjusting for age, sex, and the full kinship: ```{r} #| label: gait-assoc assoc <- association.test(geno[ok, ], Y = Y[ok], X = Xc[ok, ], method = "lmm", K = K[ok, ok], test = "score") assoc <- assoc[!is.na(assoc$p), ] ``` A **Manhattan plot** shows $-\log_{10} p$ along the genome (@fig-gwasgait): ```{r} #| label: fig-gwasgait #| fig-cap: "Manhattan plot of the mixed-model association to the XXX phenotype. The dashed line is genome-wide significance (5e-8)." #| fig-alt: "Manhattan plot with one clearly significant peak rising far above the genome-wide significance threshold on chromosome 5." #| fig-width: 8 #| fig-height: 4 manhattan(assoc, main = "XXX association (mixed model)") abline(h = -log10(5e-8), col = "#c0432f", lty = 2) ``` The scan is well-calibrated because relatedness is modelled (@fig-gait-qq): ```{r} #| label: fig-gait-qq #| fig-cap: "QQ plot of the mixed-model association p-values." #| fig-alt: "QQ plot with points on the diagonal except a short tail of genuine associations." #| fig-width: 5 #| fig-height: 5 qqplot.pvalues(assoc$p) ``` The top associated SNPs: ```{r} #| label: gait-top head(assoc[order(assoc$p), c("chr", "pos", "id", "A1", "A2", "freqA2", "p")], 10) ``` ::: {.callout-tip title="✅ Exercise 4.4 — Top hit"} 1. Locate the SNP with the smallest p-value — its `rs` id, chromosome, and position. Look it up in dbSNP. 2. Which gene lies at that locus, and what is its role in coagulation/thrombosis? 3. Compare the mixed-model result with a naïve test that ignores kinship (`association.test(geno[ok,], Y[ok], X = Xc[ok,], method = "lm")`). What changes, and why does ignoring relatedness mislead? ::: ::: {.content-hidden when-format="pdf" unless-profile="solutions"} ::: {.callout-tip collapse="true" title="Solution 4.4"} - The top hit is **rs2731672** (chr 5, ~176.8 Mb), $p \approx 2.2\times10^{-27}$, an intergenic/regulatory SNP near ***F12*** (coagulation **factor XII**, the contact-activation protease of the intrinsic pathway — directly relevant to hemostasis/thrombosis). - The modelled heritability is high ($h^2 \approx 0.82$). A naïve per-SNP linear test ignores the family correlation, so observations are not independent → p-values are anti-conservative and inflated. The mixed model uses the GRM to absorb the polygenic/family background, leaving a well-calibrated scan (the QQ plot hugs the diagonal except the genuine tail). ::: ::: ::: {.callout-important title="Ethics — genomic data privacy"} Genotype data are inherently identifying: a few dozen SNPs can re-identify an individual, and a genome reveals information about relatives who never consented. That is why the GAIT set is access-controlled and why real GWAS data live in managed-access repositories (dbGaP, EGA) under data-use agreements. Share summary statistics — never individual genotypes — and follow the study's consent and governance terms. ::: ## Further reading - Zhou & Stephens (2012), *Genetics* — the polygenic linear mixed model (`lmm.aireml`) implemented here. - [gaston vignette "linear mixed model"](https://cran.r-project.org/package=gaston) — `association.test(..., method = "lmm")`. - Soria *et al.*, the GAIT project papers, for the thrombophilia phenotype background. ::: {.content-visible when-format="html" unless-profile="solutions"} Worked solutions for the exercises are in the collapsed callouts below each exercise. ::: ::: {.content-visible when-format="html" when-profile="solutions"} Solutions to Exercises 4.1–4.4 are in the [Solutions appendix](../appendix/solutions.qmd#lab-4--gwas-gait) and in the collapsed callouts above. ::: ::: {.content-visible when-format="pdf" unless-profile="solutions"} Worked solutions for this lab are omitted from this PDF — find them as collapsed callouts in the HTML edition. ::: ## Report {#sec-report-gait} Write a laboratory report following the structure of this text. It should contain an introduction, the goals of each section, and a description of the results of each part. Include any graphics relevant to making your work understandable. **The lab should be delivered within two weeks.**