UQWS Genetic Mapping, session 1 - introduction to GWAS
Kathryn Kemper
7/6/2026
Objective
In this practical we will conduct a small GWAS in R, complete with a
basic Manhattan plot. Part 2 makes a QQ-plot by-hand, this is explained
further below.
Part 1 will conduct a small GWAS
Part 2 will plot a QQ-plot for the association study in Part
1.
You may download the data files directly to your laptops if you wish.
The data for the prac can be found in the following directory:
/data/module1/1_introToGWAS/
Use the following commands in R on the cluster to define the path to
the data. If you have downloaded the data you will need to update the
path for dir.
dir="/data/module1/1_introtoGWAS/"If you are working on the cluster you will need to save the plots
below using commands explained in the cluster information session,
i.e. if the plot is made using plot() then:
jpeg(file="myPlot.jpg") # open the jpeg.
plot(...) # call the plot.
dev.off() # close the device.
Part 1: a toy association study in R
We will now conduct a small example GWAS in R. The genotypes are from the European ancestry 1000 genomes data, and there are some simulated phenotypes.
There are three relevant files in the directory above are:
1. phenotypes: sim.phen
2. genotypes: 1kg_chr2v3_EUR.raw
3. map file: map.txt
First read in these files and investigate what they look like using
(for example) the head() and dim() functions
in R. How many individuals do you have?
# read in the genotypes, this is plink .raw format
# remove 1st 6 columns but save the IDs so we can match them with phenotype file
geno = read.table(paste0(dir,"1kg_chr2v3_EUR.raw"), header=T)
info = geno[,1:6] # save IDs
geno <- geno[,-1:-6] # remove columns from geno object
dim(geno) # there are 503 individuals, and 9,801 SNPs## [1] 503 9801 # map file
map <- read.table(paste0(dir,"map.txt"),header=T)
dim(map) # map file matches## [1] 9801 5 head(map) # information on the SNPs## chr rsID NA. position refAllele
## 1 2 rs17022634 0 100010044 C
## 2 2 rs3087385 0 100027151 C
## 3 2 rs17763718 0 100036592 T
## 4 2 rs3087395 0 100037832 C
## 5 2 rs7602535 0 100039664 C
## 6 2 rs3087399 0 100055158 C # phenotype file
pheno = read.table(paste0(dir,"sim.phen"))
dim(pheno)## [1] 503 3 sum(pheno[,1]==info[,1]) # check that the phenotypes & genotypes IDs are in the same order## [1] 503 head(pheno)## V1 V2 V3
## 1 HG00096 HG00096 -1.437370
## 2 HG00097 HG00097 -0.292070
## 3 HG00099 HG00099 -1.299560
## 4 HG00100 HG00100 -0.404354
## 5 HG00101 HG00101 6.189850
## 6 HG00102 HG00102 0.348969 hist(pheno[,3]) # check the distribution of the trait, what do you think?
A linear model for a single SNP
Let us start with a running a linear model on a single SNP. Try SNP number 491.
snp = 491
lm0 <- lm(pheno[,3] ~ geno[,snp])
summary(lm0)##
## Call:
## lm(formula = pheno[, 3] ~ geno[, snp])
##
## Residuals:
## Min 1Q Median 3Q Max
## -12.2858 -2.2242 -0.2176 2.1688 11.5730
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.5682 0.2227 -2.552 0.01101 *
## geno[, snp] 0.7544 0.2282 3.306 0.00101 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 3.47 on 501 degrees of freedom
## Multiple R-squared: 0.02135, Adjusted R-squared: 0.0194
## F-statistic: 10.93 on 1 and 501 DF, p-value: 0.001014What do the results mean? Is the SNP significantly associated with the trait?
The SNP is significantly associated with the trait, P = 0.001.
The effect of the allele on the trait is 0.7544 units and, from the model, the trait variance explained by the marker (i.e. the model \(R^2\) or ‘heritability’ of the marker) is 0.021 \(\sigma^2_P\). This is a large effect.
We can also calculate the variance explained by the marker as \(2p(1-p)a^2\), where \(p\) is the allele frequency and \(a\) is the allele effect size.
Can you see where this relationship comes from? have you seen 2p(1-p) previously?
Thus,
p = sum(geno[,snp])/(dim(geno)[1]*2)
p## [1] 0.3508946 varMarker = 2*p*(1-p)*lm0$coefficients[2]^2
h2 = varMarker / var(pheno[,3])
h2## geno[, snp]
## 0.02111313We will plot the genotypic means for the locus, to illustrate what our regression is estimating.
geno2 = geno[,snp]+rnorm(dim(geno)[1],mean=0,sd=0.05) #add a little jitter!
{ plot(pheno[,3] ~ geno2, pch=20, xlab = "genotype", ylab = "phenotype")
abline(lm0, col="light green", lty=2, lwd=2) }
A small association study
OK, so we are ready to write a simple loop in R to test each marker one-at-at-time for an association with the trait. We will use a linear model and save the corresponding p-value at each marker for the Manhattan plot.
This example also uses a for(i in array) {action} loop
in R. If you are not familiar with loops in R, try a basic loop such as
for (i in 1:10) {print(i)} to test it out. Make sure you
understand the code!
pVals <- rep(NA,ncol(geno)) # create a vector to save p-values into
for(i in 1:length(pVals)){
lm0 <- lm(pheno[,3] ~ geno[,i])
pVals[i] <- summary(lm0)$coefficients[2,4]
}
# basic Manhattan plot
map$mbp = map$pos*10^(-6)
{ plot(-log10(pVals)~map$mbp,xlab="Base pairs",pch=20)
abline (h = -log10(5*10^(-8)), col="orange", lwd=2 , lty = 2) }
The usual criteria (\(P<5x10^{-8}\)) is too conservative for this particular ‘association study’ as we didn’t really do 1M independent tests.
However, it is a simulated phenotype - we also know which loci were simulated with causal effects. Let us do another plot highlighting these SNPs, what do you make of the results?
causal = read.table(paste0(dir,"causativeLoci.txt"))
causalSNP = which(map$rsID%in%causal[,1])
{ plot(-log10(pVals)~map$mbp,xlab="Base pairs",pch=20)
abline (h = -log10(5*10^(-8)), col="orange", lwd=2 , lty = 2)
points(-log10(pVals)[causalSNP]~map$mbp[causalSNP], col="red")
abline(v=map$mbp[causalSNP],col="red") }
_________________________________________
Part 2: QQ plots
QQ plots compares two distributions by plotting their quantiles against each other. They are used in GWAS to compare the distribution of the observed p-values to their expected distribution under the null hypothesis (i.e. no association). QQ-plots are somewhat subjective but are a quick visual of the distribution and may be informative when trouble shooting. You may have seen QQ-plots previously in linear models to check the assumption of normality for the residuals.
QQ-plots involved plotting observed standardized residuals or test statistic (y − axis) on the theoretical quantiles for the data points (x-axis). The expectation when there is no association is that the data points will lie on the straight y = x line.
We are going to make a QQ-plot for our test statistics from our small GWAS by hand. There are many R packages to help with making nice QQ-plots, including those with 95% confidence intervals. I obtained some of my code for making QQ-plots using the base packages in R here.
A key piece of information here is to realise that, under the null hypothesis of no association between genotype and phenotype, the p-values follow a uniform (0,1) distribution. The code below generates quantiles of a uniform (0,1) distribution at evenly spaced probability points. So if there are 1000 SNPs, then the expected quantiles are 0.001, 0.002, 0.003, etc. The small adjustments (+0.5 and +1) are to keep the values away from the 0 and 1 boundaries.
# Calculate expectations
exp.pvalues<-(rank(pVals, ties.method="first")+.5)/(length(pVals)+1)
#Make plot
{ plot(-log10(exp.pvalues), -log10(pVals),
pch = 20, xlab="Theoretical Quartiles", ylab = "Sample Quartiles")
abline(0,1, lwd=1.3, col="red") }
# or can use a package such as 'qqman'
# install.packages("qqman")
library(qqman)## ## For example usage please run: vignette('qqman')## ## Citation appreciated but not required:## Turner, (2018). qqman: an R package for visualizing GWAS results using Q-Q and manhattan plots. Journal of Open Source Software, 3(25), 731, https://doi.org/10.21105/joss.00731.## qq(pVals)
# or one using code from Claude
library(ggplot2)
qq_with_ci <- function(pvals, ci = 0.95) {
n <- length(pvals)
# expected quantiles
expected <- -log10((rank(pvals, ties.method="first") - 0.5) / n)
observed <- -log10(sort(pvals))
# 95% CI using the beta distribution
# at rank i, p-value quantile follows Beta(i, n-i+1)
ranks <- 1:n
ci_lower <- -log10(qbeta((1 + ci)/2, ranks, n - ranks + 1))
ci_upper <- -log10(qbeta((1 - ci)/2, ranks, n - ranks + 1))
df <- data.frame(
expected = sort(expected),
observed = sort(observed),
ci_lower = sort(ci_lower),
ci_upper = sort(ci_upper)
)
ggplot(df, aes(x = expected, y = observed)) +
geom_ribbon(aes(ymin = ci_lower, ymax = ci_upper),
fill = "steelblue", alpha = 0.2) +
geom_point(size = 0.8, colour = "steelblue") +
geom_abline(slope = 1, intercept = 0,
colour = "red", linetype = "dashed", linewidth = 0.6) +
labs(x = expression(Expected ~ -log[10](p)),
y = expression(Observed ~ -log[10](p)),
title = "Q-Q plot") +
theme_minimal()
}
qq_with_ci(pVals)