3 Genetic Association: Quality Control and Stratification

In this lab we explore the steps taken in a genome-wide association (GWA) analysis: exploring a dataset, running a first association scan, discovering that it is confounded, cleaning the data, and re-running the scan. We use the maintained gaston package (Hervé Perdry & Claire Dandine-Roulland) together with base R, so the lab no longer depends on the archived GenABEL toolkit.¹

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

explore the phenotype and genotype data of a GWA dataset and summarise its quality
read minor-allele-frequency and Hardy–Weinberg equilibrium (HWE) statistics
run an initial genome-wide association scan and read a Manhattan plot
detect and remove genetic outliers using a genetic relationship matrix and PCA
apply quality-control filters and re-run the scan, and interpret why the results change

Requirements

The code in this chapter runs with gaston only. The data live in the repo under data/ge03d2.{bed,bim,fam} (genotypes) and data/ge03d2.phe (phenotypes).² The figures below are generated by executing the code.

The dataset, ge03d2, targets type 2 diabetes — a disorder characterised by high blood glucose in the context of insulin resistance and relative insulin deficiency. Unlike type 1 diabetes, the insulin resistance is generally “post-receptor”: the problem is with the cells that respond to insulin rather than with insulin production. Type 2 diabetes has a strong heritable component — having first-degree relatives with the disease substantially increases risk — and a mutation in the Islet Amyloid Polypeptide gene produces an earlier-onset, more severe form (Cho et al., 2003; Sakagashira et al., 1996; Zimmet et al., 2001). The dataset simulates some of the genetic structure of type 2 diabetes.

Lab at a glance

Goal: run and quality-control a genome-wide association scan. Prerequisites: Introduction to R; basic genetics. Time: ~2 h. Packages: gaston. Data: data/ge03d2 (PLINK bed) + data/ge03d2.phe. Deliverable: a lab report.

Start your R session, load gaston, and read the genotypes and phenotypes:

Code

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

Reading data/ge03d2.rds 
Reading data/ge03d2.bed

Code

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

geno is a bed.matrix: one row per individual, one column per marker, with annotation tables @ped (individuals) and @snps (markers). We ask gaston to compute the per-marker and per-individual summary statistics, and the HWE test:

Code

geno <- set.stats(geno, verbose = FALSE)   # call rates, MAF, heterozygosity
geno <- set.hwe(geno, verbose = FALSE)      # Hardy–Weinberg exact-test p-values
geno

A bed.matrix with 950 individuals and 7589 markers.
snps stats are set
  There are  2  SNPs with a callrate equal to zero
  There is one monomorphic SNP
ped stats are set

The individuals in @ped are in the same order as the rows of pheno, so we can use the phenotype columns directly. The autosomes are chromosomes 1–3; the X chromosome is handled separately.

Code

auto <- geno@snps$chr %in% c("1", "2", "3")

3.1 Data exploration

We first explore the features of the dataset: the number of samples, the phenotypes and their distribution, and the number and quality of the genotypes.

3.1.1 Phenotype data

The phenotype table holds the traits for each individual. str() reveals the internal structure; summary() gives per-trait distributions:

Code

str(pheno)

'data.frame':   950 obs. of  8 variables:
 $ id    : chr  "id4" "id10" "id25" "id33" ...
 $ sex   : int  0 1 0 0 0 0 1 1 0 1 ...
 $ age   : num  51.6 53.7 66 44.7 49.9 ...
 $ dm2   : int  1 1 1 1 1 1 1 1 1 1 ...
 $ height: num  152 170 165 174 157 ...
 $ weight: num  83 65.6 69.6 69.7 84.3 ...
 $ diet  : int  0 0 0 0 0 0 0 0 0 0 ...
 $ bmi   : num  36 22.6 25.6 22.9 34.1 ...

Code

summary(pheno[, c("sex", "age", "dm2", "bmi")])

      sex              age             dm2              bmi       
 Min.   :0.0000   Min.   :18.64   Min.   :0.0000   Min.   :13.49  
 1st Qu.:0.0000   1st Qu.:39.99   1st Qu.:0.0000   1st Qu.:24.28  
 Median :0.0000   Median :49.66   Median :0.0000   Median :27.94  
 Mean   :0.4895   Mean   :50.00   Mean   :0.4874   Mean   :29.68  
 3rd Qu.:1.0000   3rd Qu.:59.34   3rd Qu.:1.0000   3rd Qu.:34.11  
 Max.   :1.0000   Max.   :87.60   Max.   :1.0000   Max.   :64.22  
                                                   NAs    :2

The target trait dm2 is the binary type-2-diabetes status (0 = control, 1 = case). table() shows the case/control balance:

Code

table(dm2 = pheno$dm2)

dm2
  0   1 
487 463

✅ Exercise 3.1 — Phenotypes

What is the number of individuals in the study?
Split the traits by diabetes status (e.g. by(pheno[c("age","bmi","height","weight")], pheno$dm2, summary)) and check which traits differ most between the diabetic and normal populations.
Use a boxplot to explore graphically the relationship between the phenotypes. Which variables differ most between the diabetic and normal populations? Discuss.

Solution 3.1

The dataset has 950 individuals.
Comparing cases/controls (the by() summary and boxplots) shows bmi and, to a lesser extent, weight shift upward in the diabetic group; age is similar and height nearly identical. This is expected since BMI is a major type-2-diabetes covariate.

For example, a boxplot of body-mass index by diabetes status (Figure 3.1):

Code

boxplot(bmi ~ dm2, data = pheno,
        xlab = "Type 2 Diabetes", ylab = "Body Mass Index")

Side-by-side boxplots of body-mass index for control and type-2-diabetes groups; the diabetic group is shifted to higher BMI. — Figure 3.1: Body-mass index by type 2 diabetes status.

3.1.2 Genotype data

The @snps table now carries, for every marker, the call rate, the minor allele frequency (maf), the observed heterozygosity (hz) and the Hardy–Weinberg p-value (hwe). We summarise the autosomal markers:

Code

summary(geno@snps[auto, c("callrate", "maf", "hz", "hwe")])

    callrate           maf                hz              hwe        
 Min.   :0.0000   Min.   :0.00000   Min.   :0.0000   Min.   :0.0000  
 1st Qu.:0.9863   1st Qu.:0.05906   1st Qu.:0.1041   1st Qu.:0.2024  
 Median :0.9884   Median :0.14271   Median :0.2406   Median :0.4709  
 Mean   :0.9823   Mean   :0.18125   Mean   :0.2538   Mean   :0.4692  
 3rd Qu.:0.9905   3rd Qu.:0.28638   3rd Qu.:0.4086   3rd Qu.:0.7353  
 Max.   :0.9979   Max.   :0.50000   Max.   :0.5464   Max.   :1.0000  
                  NAs    :2         NAs    :2        NAs    :2

Minor allele frequency measures how common the rarer allele is; markers that are (nearly) monomorphic carry almost no information. Hardy–Weinberg equilibrium relates genotype frequencies to allele frequencies in a large, randomly mating population; a strong departure often flags a genotyping error — or genuine biology such as selection or population structure. Count the markers that fail typical thresholds:

Code

c(monomorphic      = sum(geno@snps$maf[auto] == 0, na.rm = TRUE),
  low_callrate     = sum(geno@snps$callrate[auto] < 0.95, na.rm = TRUE),
  hwe_p_below_1e_6 = sum(geno@snps$hwe[auto] < 1e-6, na.rm = TRUE))

     monomorphic     low_callrate hwe_p_below_1e_6 
               1               76              337

You can compare the HWE distribution between the two populations (diabetes / no diabetes) by re-computing it within each group:

Code

hwe_by <- function(rows) {
  g <- set.hwe(geno[rows, auto], verbose = FALSE)
  sum(g@snps$hwe < 1e-6, na.rm = TRUE)
}
c(controls = hwe_by(pheno$dm2 == 0),
  cases    = hwe_by(pheno$dm2 == 1))

controls    cases 
      61      360

✅ Exercise 3.2 — HWE

Do you see a difference in the number of HWE-violating markers between the two populations? Can you explain why a population might be out of HWE?

Solution 3.2

Controls: 61 markers with HWE p < 1e-6; cases: 360. Cases depart from HWE far more often.
A true association skews a marker’s genotype frequencies in cases (it is supposed to depart from HWE there), while genotyping errors and population structure cause departure in both groups. A large case-only excess is therefore a hallmark of genuine association; excess in controls points to artefacts.

3.2 Association

3.2.1 Initial scan

The target phenotype (dm2) is binary, so a first association solution is a logistic regression per marker, adjusting for age, sex and BMI. gaston::association.test() does exactly this with response = "binary"; we pass the covariates (plus an intercept column) in X. We restrict the scan to individuals with complete covariates and to the autosomes.

Code

ok  <- with(pheno, !is.na(dm2) & !is.na(age) & !is.na(sex) & !is.na(bmi))
Xc  <- cbind(1, pheno$age[ok], pheno$sex[ok], pheno$bmi[ok])
gwa1 <- association.test(geno[ok, auto], Y = pheno$dm2[ok], X = Xc,
                         method = "lm", response = "binary", test = "wald")
nrow(gwa1)

[1] 7020

The scan output has one row per marker with an association p-value. A Manhattan plot shows genome-wide significance ($-\log_{10}$ p) along the genome (Figure 3.2):

Code

manhattan(gwa1, pch = 16)

Manhattan plot of negative-log10 p-values across chromosomes 1 to 3 for the initial scan, with several peaks. — Figure 3.2: Manhattan plot of genome-wide significance (initial scan, before quality control).

The most significant markers are ranked by p-value:

Code

head(gwa1[order(gwa1$p), c("chr", "pos", "id", "p")], 10)

3.2.1.1 Is the scan well calibrated?

Before trusting any hit we check calibration with a QQ plot — observed $-\log_{10}$ p-values against those expected under the null — and the genomic inflation factor $\lambda$, the ratio of the observed median test statistic to its null expectation. $\lambda \approx 1$ means confounding is under control; $\lambda \gg 1$ signals population structure or genotyping artefacts inflating the whole scan (Figure 3.3).

Code

genomic_lambda <- function(p) {
  p <- p[is.finite(p) & p > 0]
  median(qchisq(1 - p, df = 1)) / qchisq(0.5, df = 1)
}
lambda1 <- genomic_lambda(gwa1$p)

Code

qqplot.pvalues(gwa1$p)
title(sub = sprintf("Initial scan: lambda = %.2f", lambda1))

Quantile-quantile plot of observed versus expected negative-log10 p-values for the initial association scan, with points rising above the diagonal. — Figure 3.3: QQ plot of the initial scan. Points lifting off the diagonal early indicate inflation.

The inflation factor is $\lambda = 1.22$ — above 1, a red flag that the raw scan is confounded (exactly what the quality control below will address).

Multiple testing

With thousands of markers tested at once, an uncorrected $p < 0.05$ yields hundreds of false positives. Judge hits against a genome-wide threshold — Bonferroni ($0.05/m$ for $m$ markers) or a controlled false-discovery rate (FDR, Benjamini–Hochberg via p.adjust(gwa1$p, "BH")).

✅ Exercise 3.3 — Initial scan

Can you find whether any of the markers in the table have a clinical relationship? What genome-wide significance threshold would you use for $m$ markers?

Solution 3.3

With m = 7020 markers, the Bonferroni threshold is $0.05/7020 \approx 7.1 \times 10^{-6}$; the conventional genome-wide line is $5\times10^{-8}$.
The top hit rs7903146 (a TCF7L2 variant) is the best-established T2D risk SNP (dbSNP/GWAS Catalog), so it does have a clinical relationship — a good positive control.

3.2.2 Quality control

The odd HWE distribution and probable genotyping errors are blinding our analysis. Before a second scan we clean the dataset. We drop markers with a low call rate, a very low minor allele frequency, or a gross HWE violation, and we drop individuals with a low call rate or extreme genome-wide heterozygosity (a classic signal of sample contamination or poor DNA quality).

Code

snp_ok <- geno@snps$callrate >= 0.95 &
          geno@snps$maf      >= 0.01 &
          (is.na(geno@snps$hwe) | geno@snps$hwe >= 1e-6)
c(autosomal_markers = sum(auto), pass_qc = sum(auto & snp_ok))

autosomal_markers           pass_qc 
             7020              6389

Code

het   <- geno@ped$hz
ind_ok <- geno@ped$callrate >= 0.95 &
          abs(het - mean(het)) <= 3 * sd(het)
c(individuals = nrow(geno), low_callrate = sum(geno@ped$callrate < 0.95),
  het_outliers = sum(abs(het - mean(het)) > 3 * sd(het)),
  pass_qc = sum(ind_ok))

 individuals low_callrate het_outliers      pass_qc 
         950            4            8          938

Next we look for genetic outliers — a hidden subpopulation that, mixed with the rest, inflates the whole scan. We build the genetic relationship matrix (GRM) from common autosomal markers and take its leading principal components: individuals cluster by ancestry along the first PCs.³ Figure 3.4 shows the result.

Code

gc <- geno[ok, auto & geno@snps$maf > 0.05 & geno@snps$callrate > 0.95]
standardize(gc) <- "p"
K   <- GRM(gc)
pcs <- eigen(K, symmetric = TRUE)$vectors[, 1:2]
colnames(pcs) <- c("PC1", "PC2")
km  <- kmeans(pcs, centers = 2, nstart = 50)
major <- as.integer(names(which.max(table(km$cluster))))
table(cluster = km$cluster)

cluster
  1   2 
923  25

Code

plot(pcs, col = ifelse(km$cluster == major, "grey40", "firebrick"),
     pch = 16, xlab = "PC 1", ylab = "PC 2")
legend("topright", c("main population", "outliers"),
       col = c("grey40", "firebrick"), pch = 16, bty = "n")

Scatter plot of principal component 1 versus 2 showing a large cluster and a small separated cluster of outlier individuals. — Figure 3.4: Genetic outliers in the first two principal components of the GRM. The minor cluster is the outlier subpopulation.

The minor cluster is discarded as an outlier population, as can be understood graphically from Figure 3.4.

An alternative to removing the outliers is to keep everyone and add the leading principal components as covariates in the association model (X = cbind(1, PC1, PC2, age, sex, bmi)). Both strategies control the stratification; here we remove the small outlier cluster so the effect on $\lambda$ is easy to see.

✅ Exercise 3.4 — Outliers

Re-run the k-means with a different number of centres, or colour the PCA by dm2. Does the minor cluster correspond to cases, controls, or neither?
A sex check compares the reported sex against the heterozygosity on the X chromosome (true males are hemizygous, so their X heterozygosity is ~0). Using geno@snps$chr == "X", compute per-individual X heterozygosity and look for mismatches. Why is this an important QC step?

Solution 3.4

The outlier cluster (default 2-means, 25 vs 923) is a genetic ancestry subpopulation, not the case/control split — colouring by dm2 it contains both cases and controls. So it confounds (stratifies) the scan rather than being the biological signal.
X heterozygosity flags samples where the recorded sex disagrees with the genotype (e.g. a swap or mislabelling), which would corrupt any sex-adjusted test.

We now assemble the clean dataset: the individuals in the main cluster that also pass the per-individual QC, and the markers that pass the per-marker QC.

Code

keep_ind <- which(ok)[km$cluster == major]
keep_ind <- keep_ind[ind_ok[keep_ind]]
geno_qc  <- geno[keep_ind, auto & snp_ok]
pheno_qc <- pheno[keep_ind, ]
c(individuals = nrow(geno_qc), markers = ncol(geno_qc))

individuals     markers 
        911        6389

3.2.3 Final scan

With the dataset clean — free of the stratified outliers, low-quality markers and individuals — we re-run the genome-wide scan against the type 2 diabetes phenotype (Figure 3.5):

Code

Xc2  <- cbind(1, pheno_qc$age, pheno_qc$sex, pheno_qc$bmi)
gwa2 <- association.test(geno_qc, Y = pheno_qc$dm2, X = Xc2,
                         method = "lm", response = "binary", test = "wald")
lambda2 <- genomic_lambda(gwa2$p)

Code

manhattan(gwa2, pch = 16)

Manhattan plot after quality control, with a single dominant peak and a flat background. — Figure 3.5: Manhattan plot of genome-wide significance after quality control of the dataset.

Compare these results with the first markers from the Section 3.2.1.

Code

head(gwa2[order(gwa2$p), c("chr", "pos", "id", "p")], 10)

Re-checking calibration after QC, the inflation should have fallen back towards 1 (Figure 3.6):

Code

qqplot.pvalues(gwa2$p)
title(sub = sprintf("After QC: lambda = %.2f", lambda2))

Quantile-quantile plot after quality control, with points following the diagonal except a short tail of true associations. — Figure 3.6: QQ plot after quality control — the bulk of points now track the null diagonal, with a few genuine hits departing at the tail.

Inflation dropped from $\lambda = 1.22$ (initial) to $\lambda = 1.01$ after QC: removing the stratified outliers and bad markers de-confounded the scan. The strongest surviving hit is rs7903146, a variant in the TCF7L2 gene and the best-established common risk allele for type 2 diabetes — a reassuring positive control.

✅ Exercise 3.5 — Final scan

Look up the top markers (e.g. rs7903146) in dbSNP or the GWAS Catalog. Do any have a clinical relationship to type 2 diabetes?
Can you explain, in your own words, why the results differ from the initial scan?

Solution 3.5

rs7903146 (TCF7L2) is the textbook common T2D risk variant — yes, a well documented clinical association.
The initial scan was inflated (lambda about 1.22) by the stratified ancestry outlier cluster and by poor-quality markers; after removing them the scan is recalibrated (lambda about 1.01) and the noise floor drops, so genuine hits like rs7903146 dominate instead of being masked by false positives.

3.3 Further reading

Aulchenko et al. (2007), Bioinformatics 23:1294–1296 — the original GenABEL genome-wide association framework.
The gaston vignettes — association.test, the GRM, and the mixed models used next.
PLINK2 documentation on the bed/bim/fam format.

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

3.4 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 (descriptives.marker, qtscore, check.marker, ibs); GenABEL has since been archived from CRAN. gaston reads the same data in PLINK bed format and provides the equivalent building blocks — per-marker and per-individual statistics, a genetic relationship matrix (GRM) for stratification, and a per-marker association test — which is what we use here.↩︎
ge03d2 is the simulated type-2-diabetes example that shipped with the GenABEL “getting started with GWA” tutorial (package GenABEL.data). Because GenABEL is archived, we exported it once to PLINK bed format with scripts/make-ge03d2-bed.R and committed the result, so the lab reads it with gaston and never loads GenABEL. Provenance and archived sources: GenABEL.data on the CRAN archive, the GenABEL archive, and Aulchenko et al. (2007), Bioinformatics 23(10):1294–1296.↩︎
This replaces the older Identity-By-State + multidimensional-scaling recipe. Both project a matrix of genetic similarity between every pair of individuals into a few dimensions; the eigendecomposition (PCA) of the GRM is the modern standard and is what tools such as PLINK and EIGENSTRAT use.↩︎

--- title: "Genetic Association: Quality Control and Stratification" --- In this lab we explore the steps taken in a genome-wide association (GWA) analysis: exploring a dataset, running a first association scan, discovering that it is confounded, cleaning the data, and re-running the scan. We use the maintained **gaston** package (Hervé Perdry & Claire Dandine-Roulland) together with base R, so the lab no longer depends on the archived GenABEL toolkit.[^gaston-note] [^gaston-note]: Earlier editions of this lab used **GenABEL** (`descriptives.marker`, `qtscore`, `check.marker`, `ibs`); GenABEL has since been archived from CRAN. `gaston` reads the same data in PLINK **bed** format and provides the equivalent building blocks — per-marker and per-individual statistics, a genetic relationship matrix (GRM) for stratification, and a per-marker association test — which is what we use here. By the end of this laboratory you will be able to: - explore the phenotype and genotype data of a GWA dataset and summarise its quality - read minor-allele-frequency and Hardy–Weinberg equilibrium (HWE) statistics - run an initial genome-wide association scan and read a Manhattan plot - detect and remove genetic outliers using a genetic relationship matrix and PCA - apply quality-control filters and re-run the scan, and interpret why the results change ::: {.callout-note title="Requirements"} The code in this chapter runs with **gaston** only. The data live in the repo under `data/ge03d2.{bed,bim,fam}` (genotypes) and `data/ge03d2.phe` (phenotypes).[^data-note] The figures below are generated by executing the code. ::: [^data-note]: `ge03d2` is the simulated type-2-diabetes example that shipped with the GenABEL "getting started with GWA" tutorial (package `GenABEL.data`). Because GenABEL is archived, we exported it once to PLINK bed format with `scripts/make-ge03d2-bed.R` and committed the result, so the lab reads it with `gaston` and never loads GenABEL. Provenance and archived sources: [GenABEL.data on the CRAN archive](https://cran.r-project.org/src/contrib/Archive/GenABEL.data/), [the GenABEL archive](https://cran.r-project.org/src/contrib/Archive/GenABEL/), and Aulchenko *et al.* (2007), *Bioinformatics* **23**(10):1294–1296. The dataset, `ge03d2`, targets **type 2 diabetes** — a disorder characterised by high blood glucose in the context of insulin resistance and relative insulin deficiency. Unlike type 1 diabetes, the insulin resistance is generally "post-receptor": the problem is with the cells that respond to insulin rather than with insulin production. Type 2 diabetes has a strong heritable component — having first-degree relatives with the disease substantially increases risk — and a mutation in the *Islet Amyloid Polypeptide* gene produces an earlier-onset, more severe form [@Zimmet2001; @Sakagashira1996; @Cho2003]. The dataset simulates some of the genetic structure of type 2 diabetes. ```{r} #| label: assoc-setup #| include: false source("R/setup.R") options(width = 75) ``` ::: {.callout-note appearance="minimal" title="Lab at a glance"} **Goal:** run and quality-control a genome-wide association scan. **Prerequisites:** [Introduction to R](intro-r.qmd); basic genetics. **Time:** ~2 h. **Packages:** `gaston`. **Data:** `data/ge03d2` (PLINK bed) + `data/ge03d2.phe`. **Deliverable:** a lab report. ::: Start your R session, load `gaston`, and read the genotypes and phenotypes: ```{r} #| label: load-data #| message: false library(gaston) geno <- read.bed.matrix("data/ge03d2") pheno <- read.table("data/ge03d2.phe", header = TRUE) ``` `geno` is a **bed.matrix**: one row per individual, one column per marker, with annotation tables `@ped` (individuals) and `@snps` (markers). We ask `gaston` to compute the per-marker and per-individual summary statistics, and the HWE test: ```{r} #| label: set-stats geno <- set.stats(geno, verbose = FALSE) # call rates, MAF, heterozygosity geno <- set.hwe(geno, verbose = FALSE) # Hardy–Weinberg exact-test p-values geno ``` The individuals in `@ped` are in the same order as the rows of `pheno`, so we can use the phenotype columns directly. The autosomes are chromosomes 1–3; the X chromosome is handled separately. ```{r} #| label: autosomes auto <- geno@snps$chr %in% c("1", "2", "3") ``` ## Data exploration {#sec-data-exploration} We first explore the features of the dataset: the number of samples, the phenotypes and their distribution, and the number and quality of the genotypes. ### Phenotype data {#sec-phenotype-data} The phenotype table holds the traits for each individual. `str()` reveals the internal structure; `summary()` gives per-trait distributions: ```{r} #| label: pheno-summary str(pheno) summary(pheno[, c("sex", "age", "dm2", "bmi")]) ``` The target trait `dm2` is the binary type-2-diabetes status (0 = control, 1 = case). `table()` shows the case/control balance: ```{r} #| label: dm2-table table(dm2 = pheno$dm2) ``` ::: {.callout-tip title="✅ Exercise 3.1 — Phenotypes"} 1. What is the number of individuals in the study? 2. Split the traits by diabetes status (e.g. `by(pheno[c("age","bmi","height","weight")], pheno$dm2, summary)`) and check which traits differ most between the diabetic and normal populations. 3. Use a boxplot to explore graphically the relationship between the phenotypes. Which variables differ most between the diabetic and normal populations? Discuss. ::: ::: {.content-hidden when-format="pdf" unless-profile="solutions"} ::: {.callout-tip collapse="true" title="Solution 3.1"} - The dataset has **950 individuals**. - Comparing cases/controls (the `by()` summary and boxplots) shows **bmi** and, to a lesser extent, **weight** shift upward in the diabetic group; **age** is similar and **height** nearly identical. This is expected since BMI is a major type-2-diabetes covariate. ::: ::: For example, a boxplot of body-mass index by diabetes status (@fig-bmi-box): ```{r} #| label: fig-bmi-box #| fig-cap: "Body-mass index by type 2 diabetes status." #| fig-alt: "Side-by-side boxplots of body-mass index for control and type-2-diabetes groups; the diabetic group is shifted to higher BMI." #| fig-width: 5 #| fig-height: 4 boxplot(bmi ~ dm2, data = pheno, xlab = "Type 2 Diabetes", ylab = "Body Mass Index") ``` ### Genotype data {#sec-genotype-data} The `@snps` table now carries, for every marker, the call rate, the minor allele frequency (`maf`), the observed heterozygosity (`hz`) and the Hardy–Weinberg p-value (`hwe`). We summarise the autosomal markers: ```{r} #| label: marker-summary summary(geno@snps[auto, c("callrate", "maf", "hz", "hwe")]) ``` **Minor allele frequency** measures how common the rarer allele is; markers that are (nearly) monomorphic carry almost no information. **Hardy–Weinberg equilibrium** relates genotype frequencies to allele frequencies in a large, randomly mating population; a strong departure often flags a genotyping error — or genuine biology such as selection or population structure. Count the markers that fail typical thresholds: ```{r} #| label: marker-counts c(monomorphic = sum(geno@snps$maf[auto] == 0, na.rm = TRUE), low_callrate = sum(geno@snps$callrate[auto] < 0.95, na.rm = TRUE), hwe_p_below_1e_6 = sum(geno@snps$hwe[auto] < 1e-6, na.rm = TRUE)) ``` You can compare the HWE distribution between the two populations (diabetes / no diabetes) by re-computing it within each group: ```{r} #| label: hwe-two hwe_by <- function(rows) { g <- set.hwe(geno[rows, auto], verbose = FALSE) sum(g@snps$hwe < 1e-6, na.rm = TRUE) } c(controls = hwe_by(pheno$dm2 == 0), cases = hwe_by(pheno$dm2 == 1)) ``` ::: {.callout-tip title="✅ Exercise 3.2 — HWE"} Do you see a difference in the number of HWE-violating markers between the two populations? Can you explain why a population might be out of HWE? ::: ::: {.content-hidden when-format="pdf" unless-profile="solutions"} ::: {.callout-tip collapse="true" title="Solution 3.2"} - Controls: **61** markers with HWE p < 1e-6; cases: **360**. Cases depart from HWE far more often. - A true association skews a marker's genotype frequencies in cases (it is *supposed* to depart from HWE there), while genotyping errors and population structure cause departure in both groups. A large case-only excess is therefore a hallmark of genuine association; excess in controls points to artefacts. ::: ::: ## Association {#sec-association} ### Initial scan {#sec-initial-scan} The target phenotype (`dm2`) is binary, so a first association solution is a **logistic regression per marker**, adjusting for age, sex and BMI. `gaston::association.test()` does exactly this with `response = "binary"`; we pass the covariates (plus an intercept column) in `X`. We restrict the scan to individuals with complete covariates and to the autosomes. ```{r} #| label: scan-initial ok <- with(pheno, !is.na(dm2) & !is.na(age) & !is.na(sex) & !is.na(bmi)) Xc <- cbind(1, pheno$age[ok], pheno$sex[ok], pheno$bmi[ok]) gwa1 <- association.test(geno[ok, auto], Y = pheno$dm2[ok], X = Xc, method = "lm", response = "binary", test = "wald") nrow(gwa1) ``` The scan output has one row per marker with an association p-value. A **Manhattan plot** shows genome-wide significance ($-\log_{10}$ p) along the genome (@fig-manhattan-initial): ```{r} #| label: fig-manhattan-initial #| fig-cap: "Manhattan plot of genome-wide significance (initial scan, before quality control)." #| fig-alt: "Manhattan plot of negative-log10 p-values across chromosomes 1 to 3 for the initial scan, with several peaks." #| fig-width: 6 #| fig-height: 4 manhattan(gwa1, pch = 16) ``` The most significant markers are ranked by p-value: ```{r} #| label: scan-results-gwa1 head(gwa1[order(gwa1$p), c("chr", "pos", "id", "p")], 10) ``` #### Is the scan well calibrated? {#sec-qq-initial} Before trusting any hit we check **calibration** with a QQ plot — observed $-\log_{10}$ p-values against those expected under the null — and the **genomic inflation factor** $\lambda$, the ratio of the observed median test statistic to its null expectation. $\lambda \approx 1$ means confounding is under control; $\lambda \gg 1$ signals population structure or genotyping artefacts inflating the whole scan (@fig-qq-initial). ```{r} #| label: lambda-fn genomic_lambda <- function(p) { p <- p[is.finite(p) & p > 0] median(qchisq(1 - p, df = 1)) / qchisq(0.5, df = 1) } lambda1 <- genomic_lambda(gwa1$p) ``` ```{r} #| label: fig-qq-initial #| fig-cap: "QQ plot of the initial scan. Points lifting off the diagonal early indicate inflation." #| fig-alt: "Quantile-quantile plot of observed versus expected negative-log10 p-values for the initial association scan, with points rising above the diagonal." #| fig-width: 5 #| fig-height: 5 qqplot.pvalues(gwa1$p) title(sub = sprintf("Initial scan: lambda = %.2f", lambda1)) ``` The inflation factor is $\lambda = `r round(lambda1, 2)`$ — above 1, a red flag that the raw scan is confounded (exactly what the quality control below will address). ::: {.callout-important title="Multiple testing"} With thousands of markers tested at once, an uncorrected $p < 0.05$ yields hundreds of false positives. Judge hits against a genome-wide threshold — Bonferroni ($0.05/m$ for $m$ markers) or a controlled false-discovery rate (FDR, Benjamini–Hochberg via `p.adjust(gwa1$p, "BH")`). ::: ::: {.callout-tip title="✅ Exercise 3.3 — Initial scan"} Can you find whether any of the markers in the table have a clinical relationship? What genome-wide significance threshold would you use for $m$ markers? ::: ::: {.content-hidden when-format="pdf" unless-profile="solutions"} ::: {.callout-tip collapse="true" title="Solution 3.3"} - With **m = 7020** markers, the Bonferroni threshold is $0.05/7020 \approx 7.1 \times 10^{-6}$; the conventional genome-wide line is $5\times10^{-8}$. - The top hit **rs7903146** (a *TCF7L2* variant) is the best-established T2D risk SNP (dbSNP/GWAS Catalog), so it does have a clinical relationship — a good positive control. ::: ::: ### Quality control {#sec-quality-control} The odd HWE distribution and probable genotyping errors are blinding our analysis. Before a second scan we clean the dataset. We drop **markers** with a low call rate, a very low minor allele frequency, or a gross HWE violation, and we drop **individuals** with a low call rate or extreme genome-wide heterozygosity (a classic signal of sample contamination or poor DNA quality). ```{r} #| label: qc-markers snp_ok <- geno@snps$callrate >= 0.95 & geno@snps$maf >= 0.01 & (is.na(geno@snps$hwe) | geno@snps$hwe >= 1e-6) c(autosomal_markers = sum(auto), pass_qc = sum(auto & snp_ok)) ``` ```{r} #| label: qc-individuals het <- geno@ped$hz ind_ok <- geno@ped$callrate >= 0.95 & abs(het - mean(het)) <= 3 * sd(het) c(individuals = nrow(geno), low_callrate = sum(geno@ped$callrate < 0.95), het_outliers = sum(abs(het - mean(het)) > 3 * sd(het)), pass_qc = sum(ind_ok)) ``` Next we look for **genetic outliers** — a hidden subpopulation that, mixed with the rest, inflates the whole scan. We build the **genetic relationship matrix** (GRM) from common autosomal markers and take its leading **principal components**: individuals cluster by ancestry along the first PCs.[^pca-note] @fig-pca shows the result. [^pca-note]: This replaces the older Identity-By-State + multidimensional-scaling recipe. Both project a matrix of genetic similarity between every pair of individuals into a few dimensions; the eigendecomposition (PCA) of the GRM is the modern standard and is what tools such as PLINK and EIGENSTRAT use. ```{r} #| label: grm-pca gc <- geno[ok, auto & geno@snps$maf > 0.05 & geno@snps$callrate > 0.95] standardize(gc) <- "p" K <- GRM(gc) pcs <- eigen(K, symmetric = TRUE)$vectors[, 1:2] colnames(pcs) <- c("PC1", "PC2") km <- kmeans(pcs, centers = 2, nstart = 50) major <- as.integer(names(which.max(table(km$cluster)))) table(cluster = km$cluster) ``` ```{r} #| label: fig-pca #| fig-cap: "Genetic outliers in the first two principal components of the GRM. The minor cluster is the outlier subpopulation." #| fig-alt: "Scatter plot of principal component 1 versus 2 showing a large cluster and a small separated cluster of outlier individuals." #| fig-width: 5 #| fig-height: 5 plot(pcs, col = ifelse(km$cluster == major, "grey40", "firebrick"), pch = 16, xlab = "PC 1", ylab = "PC 2") legend("topright", c("main population", "outliers"), col = c("grey40", "firebrick"), pch = 16, bty = "n") ``` The minor cluster is discarded as an outlier population, as can be understood graphically from @fig-pca. ::: {.callout-note appearance="simple"} An alternative to *removing* the outliers is to keep everyone and add the leading principal components as covariates in the association model (`X = cbind(1, PC1, PC2, age, sex, bmi)`). Both strategies control the stratification; here we remove the small outlier cluster so the effect on $\lambda$ is easy to see. ::: ::: {.callout-tip title="✅ Exercise 3.4 — Outliers"} 1. Re-run the k-means with a different number of centres, or colour the PCA by `dm2`. Does the minor cluster correspond to cases, controls, or neither? 2. A **sex check** compares the reported sex against the heterozygosity on the X chromosome (true males are hemizygous, so their X heterozygosity is ~0). Using `geno@snps$chr == "X"`, compute per-individual X heterozygosity and look for mismatches. Why is this an important QC step? ::: ::: {.content-hidden when-format="pdf" unless-profile="solutions"} ::: {.callout-tip collapse="true" title="Solution 3.4"} - The outlier cluster (default 2-means, **25 vs 923**) is a *genetic ancestry* subpopulation, not the case/control split — colouring by `dm2` it contains both cases and controls. So it confounds (stratifies) the scan rather than being the biological signal. - X heterozygosity flags samples where the recorded sex disagrees with the genotype (e.g. a swap or mislabelling), which would corrupt any sex-adjusted test. ::: ::: We now assemble the clean dataset: the individuals in the main cluster that also pass the per-individual QC, and the markers that pass the per-marker QC. ```{r} #| label: clean-dataset keep_ind <- which(ok)[km$cluster == major] keep_ind <- keep_ind[ind_ok[keep_ind]] geno_qc <- geno[keep_ind, auto & snp_ok] pheno_qc <- pheno[keep_ind, ] c(individuals = nrow(geno_qc), markers = ncol(geno_qc)) ``` ### Final scan {#sec-final-scan} With the dataset clean — free of the stratified outliers, low-quality markers and individuals — we re-run the genome-wide scan against the type 2 diabetes phenotype (@fig-manhattan-qc): ```{r} #| label: scan-final Xc2 <- cbind(1, pheno_qc$age, pheno_qc$sex, pheno_qc$bmi) gwa2 <- association.test(geno_qc, Y = pheno_qc$dm2, X = Xc2, method = "lm", response = "binary", test = "wald") lambda2 <- genomic_lambda(gwa2$p) ``` ```{r} #| label: fig-manhattan-qc #| fig-cap: "Manhattan plot of genome-wide significance after quality control of the dataset." #| fig-alt: "Manhattan plot after quality control, with a single dominant peak and a flat background." #| fig-width: 6 #| fig-height: 4 manhattan(gwa2, pch = 16) ``` Compare these results with the first markers from the @sec-initial-scan. ```{r} #| label: tablegwa2 head(gwa2[order(gwa2$p), c("chr", "pos", "id", "p")], 10) ``` Re-checking calibration after QC, the inflation should have fallen back towards 1 (@fig-qq-qc): ```{r} #| label: fig-qq-qc #| fig-cap: "QQ plot after quality control — the bulk of points now track the null diagonal, with a few genuine hits departing at the tail." #| fig-alt: "Quantile-quantile plot after quality control, with points following the diagonal except a short tail of true associations." #| fig-width: 5 #| fig-height: 5 qqplot.pvalues(gwa2$p) title(sub = sprintf("After QC: lambda = %.2f", lambda2)) ``` Inflation dropped from $\lambda = `r round(lambda1, 2)`$ (initial) to $\lambda = `r round(lambda2, 2)`$ after QC: removing the stratified outliers and bad markers de-confounded the scan. The strongest surviving hit is **rs7903146**, a variant in the *TCF7L2* gene and the best-established common risk allele for type 2 diabetes — a reassuring positive control. ::: {.callout-tip title="✅ Exercise 3.5 — Final scan"} 1. Look up the top markers (e.g. rs7903146) in [dbSNP](https://www.ncbi.nlm.nih.gov/snp/) or the GWAS Catalog. Do any have a clinical relationship to type 2 diabetes? 2. Can you explain, in your own words, why the results differ from the initial scan? ::: ::: {.content-hidden when-format="pdf" unless-profile="solutions"} ::: {.callout-tip collapse="true" title="Solution 3.5"} - **rs7903146 (TCF7L2)** is the textbook common T2D risk variant — yes, a well documented clinical association. - The initial scan was inflated (lambda about 1.22) by the stratified ancestry outlier cluster and by poor-quality markers; after removing them the scan is recalibrated (lambda about 1.01) and the noise floor drops, so genuine hits like rs7903146 dominate instead of being masked by false positives. ::: ::: ## Further reading - Aulchenko *et al.* (2007), *Bioinformatics* **23**:1294–1296 — the original GenABEL genome-wide association framework. - The [gaston vignettes](https://cran.r-project.org/package=gaston) — `association.test`, the GRM, and the mixed models used next. - [PLINK2 documentation](https://www.cog-genomics.org/plink/2.0/) on the bed/bim/fam format. ::: {.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 3.1–3.5 are in the [Solutions appendix](../appendix/solutions.qmd#lab-3--genetic-association) 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-assoc} 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.**