library(ggplot2)
library(reshape2)

Objective

In this practical we will take a look at some genotype files from the 1000 genomes project. You can choose part 1, 2, 3 or 4. If this is all new, start with Part 1.

Part 1 will investigate \(r^2\) as a measure of LD between SNPs and create a LD-matrix plot using LDlink.
Part 2 will investigate the power to detect loci.
Part 3 will construct a small toy example of a genomic-relationship matrix (GRM) and also look at some properties of a GRM computed with GCTA.
Part 4 will perform a simple PCA analysis in R

all the data used in this practical is publically available, you may download the data if you wish.

The data for the prac can be found in the following directory:
/data/module1/2_studyDesignPrac/

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/2_studyDesignPrac/"

Again, if you are working on the cluster you will need to save the plots below using the same commands in the first prac, i.e. if the plot is made using plot() then:

jpeg(file="myPlot.jpg",type='cairo') # open the jpeg.

plot(...) # call the plot.

dev.off() # close the device.

OR, if the plot is made using ggplot() then:

savePlot = ggplot() # save the plot as an object.

ggsave(savePlot,file="myPlot.jpg") # use ggsave() to save the object.

Part 1: LD between loci

Data used in this part 1 of the practical are:

1kg_chr2v3_EUR.raw → PLINK (text) file for 1000G european ancestry individuals.

1kg_chr2v3_AFR2.raw → PLINK (text) file for 1000G african ancestry individuals.

map.txt → corresponding map for the above data; chromosome, marker name (rsID) and position in base pairs.

LD (linkage disequilbrium) is the non-random association between genetic markers. It is created by broken down by recombination, and influenced by population demographic events. In this first part of the practical, we will use European and African ancestry genomes from the 1000 genomes project to investigate LD in a small region on chromosome 2.

The two genotype files are 1kg_chr2v3_EUR.raw (Europeans) and 1kg_chr2v3_AFR2.raw (Africans). I have converted the binary PLINK format (which cannot easily be read into R or interpreted without using software) into a text file format with the --recode A option so that the alleles are coded as 0, 1 or 2 copies of the reference allele.

Read in the files and take a look at the format.

  # read in the EUR genotypes & take a look
  EUR = read.table(paste0(dir,"1kg_chr2v3_EUR.raw"), header=T)
  head(EUR[,1:10])

##       FID     IID PAT MAT SEX PHENOTYPE rs17022634_C rs3087385_C rs17763718_T
## 1 HG00096 HG00096   0   0   0        -9            0           0            0
## 2 HG00097 HG00097   0   0   0        -9            1           0            0
## 3 HG00099 HG00099   0   0   0        -9            0           0            1
## 4 HG00100 HG00100   0   0   0        -9            0           0            0
## 5 HG00101 HG00101   0   0   0        -9            1           0            0
## 6 HG00102 HG00102   0   0   0        -9            1           0            0
##   rs3087395_C
## 1           0
## 2           0
## 3           0
## 4           0
## 5           0
## 6           0

  # we don't need information on individuals so we'll trim these columns
  EUR <- EUR[,-1:-6]
  head(EUR[,1:10])

##   rs17022634_C rs3087385_C rs17763718_T rs3087395_C rs7602535_C rs3087399_C
## 1            0           0            0           0           0           0
## 2            1           0            0           0           0           0
## 3            0           0            1           0           0           0
## 4            0           0            0           0           0           0
## 5            1           0            0           0           0           0
## 6            1           0            0           0           0           1
##   rs3087403_T rs28369955_G rs28369946_G rs10191001_C
## 1           1            0            0            0
## 2           1            0            0            0
## 3           0            0            0            0
## 4           2            0            0            0
## 5           1            0            0            0
## 6           0            0            0            0

  dim(EUR)            # 503 rows = number of individuals, 9801 columns = number of SNPs

## [1]  503 9801

  map  <- read.table(paste0(dir,"map.txt"),header=T)
  dim(map)            # 9801 rows = number of SNPs, matches above yay!

## [1] 9801    5

  AFR = read.table(paste0(dir,"1kg_chr2v3_AFR2.raw"),header=T)
  AFR <- AFR[,-1:-6]  # trim again
  dim(AFR)            # there are 660 individuals, and 9,801 SNPs

## [1]  660 9801

Calculating allele frequencies

Now we’ll check the allele frequencies in the samples. Make sure you understand the code. If necessary run subsections of the code and use head() or dim() to understand what the functions are doing.

   # first calculate the frequency at a single SNP, e.g. the first one
   snp = EUR[,1]
   snp                  # check the format, coded 0, 1 and 2

##   [1] 0 1 0 0 1 1 1 1 0 0 0 0 0 0 1 0 0 1 1 0 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0
##  [38] 0 0 0 0 0 1 1 0 0 0 0 1 0 0 0 0 0 1 0 0 1 0 0 1 0 0 0 0 0 0 0 0 1 0 0 1 2
##  [75] 0 1 0 0 0 1 0 1 0 0 0 1 0 0 1 1 0 0 0 0 1 0 0 0 0 0 0 0 1 1 1 1 0 0 0 0 1
## [112] 1 0 0 0 0 0 0 0 1 0 1 1 0 0 1 1 0 0 0 0 0 0 0 0 0 0 1 1 0 1 0 0 1 0 1 0 0
## [149] 1 1 0 0 1 0 0 0 1 0 0 0 1 0 1 0 1 0 0 0 1 0 1 1 0 0 0 0 0 0 1 0 0 1 0 0 1
## [186] 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1
## [223] 0 0 1 0 1 0 1 0 0 0 1 1 0 1 0 0 1 0 1 0 0 0 1 0 0 1 0 1 0 1 0 1 1 1 0 0 1
## [260] 0 1 0 1 1 0 1 2 0 1 0 0 0 0 0 0 0 0 1 0 0 0 1 0 0 0 1 0 0 0 0 0 0 1 1 0 0
## [297] 0 0 0 0 0 1 0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 1 1 0 0 0 1 0 0 0 0 1 0 0
## [334] 1 0 0 0 0 0 0 0 1 1 1 0 0 0 0 0 0 0 0 0 1 0 0 2 0 0 1 1 1 0 0 0 0 0 0 1 0
## [371] 0 0 0 0 0 1 0 1 0 0 0 1 0 1 0 0 0 1 1 0 0 0 0 0 0 0 1 0 0 0 1 0 1 0 1 1 0
## [408] 0 0 0 1 0 0 0 1 1 1 0 0 2 0 1 0 0 0 0 0 0 0 1 0 1 0 0 1 0 0 0 0 1 0 0 0 0
## [445] 0 1 1 1 0 0 0 1 0 1 1 0 0 0 0 0 2 0 0 0 0 0 0 0 1 0 0 0 0 0 1 0 1 0 0 0 0
## [482] 0 0 0 0 0 1 0 0 0 0 0 0 0 1 0 0 2 1 0 0 0 0

   freq = mean(snp)/2   # why divide by 2?
   freq

## [1] 0.1431412

   # now all SNPs
   EURp = colMeans(EUR)/2
   AFRp = colMeans(AFR)/2
   plot(AFRp ~ EURp, pch = 20, col="dark blue", 
      xlab = "allele frequency, EUR", ylab = "allele frequency, AFR")

   # we can also see the data has been filtered for MAF 0.01 in EUR, but not in AFR
   range(EURp)

## [1] 0.01093439 0.50000000

   range(AFRp)

## [1] 0.0000000 0.9810606

   mafSNP = AFRp>0.01
   
   # add the frequencies as columns to the map file - for convenience
   map = cbind(map, EURp, AFRp)

LD between loci

We are going to calculate LD for a marker with all surrounding markers, either in the EUR or AFR ancestry groups. I have chosen rs77732099, which is a marker at about 135 Mbp. I will calculate the \(r^2\) or squared correlation with the marker rs77732099 and all other markers within a 4Mbp region.

   marker = "rs77732099"
   SNP = which(map$rsID==marker)
   map[SNP,]      # check the SNP is not monomorphic

##              chr       rsID NA.  position refEUR      EURp       AFRp
## rs77732099_G   2 rs77732099   0 135866750      G 0.1023857 0.05227273

   # identify SNPs in 4Mbp region
   start = map[SNP,"position"] - 2*10^6
   stop = map[SNP,"position"] + 2*10^6
   region = map$position>start & map$position < stop
   testSNPs = which (region & mafSNP)
   
   # write a small loop to calculate LD with all SNP in the region
   LD_EUR = NULL ; LD_AFR = NULL
   for (i in testSNPs) {
    LD_EUR = c(LD_EUR,cor(EUR[,SNP],EUR[,i])^2)
    LD_AFR = c(LD_AFR,cor(AFR[,SNP],AFR[,i])^2)
   }
   
   # plot the result
   map$mbp = map$pos/10^6
   { plot(LD_EUR ~ map$mbp[testSNPs], pch = 20, col="orange", 
          xlab = "linkage disequilbrium, chr2 Mbp",
          ylim = range(c(LD_EUR,LD_AFR)))
     points(LD_AFR ~ map$mbp[testSNPs], pch=20, col="dark blue") 
     legend("topright", legend = c("EUR","AFR"), fill = c("orange","dark blue")) 
     abline(v=map$mbp[SNP], col="red", lty = 2) }

In practice, it is very easy in PLINK to calculate allele frequencies and/or LD for genotype files in binary format. We will look at doing this in the genotype cleaning prac. However, its useful to see how to do it manually - sometimes this is needed when trouble shooting or checking the output from PLINK.

An LD-matrix plot

We can see extensive LD in this region. The region was not a random choice, it is near the LCT (or lactase) gene, which confers lactase persistence in humans. Lactase persistence is where lactose (the primary carbohydrate found in cows milk) can be digested into adulthood. It has been under selection. A single SNP (rs4988235) confers lactase persistence in European-descent populations, with it being one of several SNPs conferring lactase persistence in African ancestries.

An LD-matrix is a square or triangle with colour coded pair-wise LD metrics. If the matrix is square, \(D'\) will be on one half of the matrix and \(r^2\) on the other. We will now use a web-based tool LDlink to investigate LD in this region using an LD-matrix plot.
1. From the home page, click on the ‘documentation’ page and then scroll down to the ‘LDmatrix’ module for the requirements.
2. The input is a list of up to 300 rsID numbers. Use R to select up to 300 SNP to input into the tool.
3. Go back to the homepage and navigate to the LDmatrix tool and create your LD-matrix plot. It should only take a few moments to produce the plot.

The LDlink tool has many modules and tools for GWAS which may be useful. There is also an R package called LDlinkR.

   length(testSNPs) # too many, select about 300 in the middle

## [1] 661

   testSNPs = testSNPs[201:500]
   write.table(file="rsIDs.txt",map[testSNPs,"rsID"],quote=FALSE,row.names=FALSE,col.names=FALSE)

In case that your plot doesn’t work in real time, here is one I generated using all the EUR populations (CEU, TSI, FIN, GBR, IBS).

An LD-matrix near LCT in European-ancestry samples

Part 2: Power to detect loci

This part of the practical doesn’t use any data, you may run it directly on your laptop or run it on the server and download any plots.

Knowledge of the factors affecting statistical power to detect loci is essential for understanding and interpreting GWAS results. We will use and expand upon the simple power calculations for quantitative traits found here. Please read.

As outlined in the wiki, the factors affecting power are:
1. \(N\), the sample size
2. \(h^2\), the ‘heritability’ of the variant. This is the proportion of phenotypic variance explained by the locus. We will expand on this parameter to include allele frequency and the effect size, as well as LD between causal variants and genotyped SNP below.
3. \(\alpha\), the significance level

First, let us try a simple example. What is the power to detect a loci explaining 1% of the phenotypic variance (i.e. heritability of the marker is 0.01) with a sample size of 500 at genome-wide significance (\(\alpha = 5x10^{-8}\))?

    N = 500                         # sample size
    alpha = 0.05/(1*10^6)           # significance level
    h2 = 0.01                       # heritability of the locus
    threshold = qchisq(alpha, df = 1, lower.tail = FALSE)
    power = pchisq(threshold, df = 1, lower.tail = FALSE, ncp = N * h2)
    power

## [1] 0.0006516719

Ouch, pretty poor! Less than 1% power to detect a loci with this sample size at \(P<5x10^{-8}\). How big does the sample size need to be to detect such as locus with 80% power?

    N = c(500,seq(1000,5000,500)) #guess some sample sizes
    power = pchisq(threshold, df = 1, lower.tail = FALSE, ncp = N * h2)
    power

##  [1] 0.0006516719 0.0110387240 0.0572452426 0.1637468872 0.3258829115
##  [6] 0.5103374868 0.6789516787 0.8087351921 0.8956038978 0.9473577969

    { plot(power ~ N, type="b", main = "sample size to detect loci explaining 1% of phenotypic variance")
      abline(h=0.8, col="orange", lty = 2, lwd = 2) }

So that’s about 4,000 (unrelated) samples with genotypes and phenotypes to detect a genotyped locus explaining 1% of the phenotypic variance. Now, what happens if you haven’t genotyped the ‘causal’ locus and you need to consider LD between the locus and the genotyped SNP?

Given a sample size of 4,000; how is the power to detect this locus reduced when the LD between the causal variant and the genotyped SNP decreases?

   N = 4000
   r2 = seq(0.05,1,0.05)
   h2 = 0.01*r2
   power = pchisq(threshold, df = 1, lower.tail = FALSE, ncp = N * h2)
   
   { plot(power ~ r2, type="b", main = "reduction in power due to incomplete LD with causal variant")
      abline(h=0.8, col="orange", lty = 2, lwd = 2)}

The variance explained by a locus is dependent on allele frequency, and it can be calculated as \(2p(1-p)a^2\), where \(p\) is the allele frequency and \(a\) is the effect size in trait units.

Use the power calculator to show the reduction in power by changing the allele frequency of an allele when the allele explains 0.01 of \(\sigma^2_P\) at a frequency of \(p = 0.5\)

   N = 4000
   a = sqrt(0.01/(2*0.5^2))   # h2SNP = 2p(1-p)a^2; backsolve for the effect size a
   p = seq(0.01,0.99,0.01)    # range of allele frequencies
   h2 = 2*p*(1-p)*a^2         # variance explained by the locus
   power = pchisq(threshold, df = 1, lower.tail = FALSE, ncp = N * h2)
   
   plot(power ~ p, type="b", main = "influence of allele frequency to detect loci", las = 1)

Part 3: A simple relationship matrix

Data used in the part of the practical are:

1kg_rand_EUR.raw. → same file as Part1, european ancestry genotypes from 1000G.

families.grm.gz. → a GCTA genomic relationship matrix (GRM) file in g-zipped format.

familes.grm.id. → a GCTA file of sample IDs corresponding to the GRM above.

The concept of a genomic relationship matrix (GRM) is useful when dealing with relatives and more complex GWAS models. It is a square matrix that, in it’s most basic form, is calculated by:
\[ GRM = WW'/n \] where \(W\) is a matrix of scaled genotypes (scaled so that each SNP has mean 0 and variance 1), and \(n\) is the number of markers.

making the GRM, toy example

We will make a simple GRM and plot it in unrelated individuals. First read in the data, and check MAF > 0.01. why is it important to ensure that MAF > 0.01?

  # read in the EUR genotypes
  EUR2 = read.table(paste0(dir,"1kg_rand_EUR.raw"), header=T)
  dim(EUR2)           # 503 individuals, 457 SNP

## [1] 503 457

  head(EUR2[,1:10])   # same format as above

##       FID     IID PAT MAT SEX PHENOTYPE rs4908438_G rs10737914_T rs6704445_A
## 1 HG00096 HG00096   0   0   0        -9           0            1           1
## 2 HG00097 HG00097   0   0   0        -9           0            2           1
## 3 HG00099 HG00099   0   0   0        -9           0            0           1
## 4 HG00100 HG00100   0   0   0        -9           0            1           0
## 5 HG00101 HG00101   0   0   0        -9           0            2           1
## 6 HG00102 HG00102   0   0   0        -9           0            1           1
##   rs34824071_A
## 1            1
## 2            1
## 3            1
## 4            0
## 5            0
## 6            2

  EUR2 = EUR2[,-1:-6] # remove ids
  EUR2p = colMeans(EUR2)/2
  range(EUR2p)        # filtered for maf

## [1] 0.01093439 0.49900596

Making the GRM is simple on this small scale. However, as we have only used very few SNPs (n = 451), the properties of the GRM are more variable than when using large numbers of variants. Constructing the GRM can be very time consuming for large numbers of individuals with large numbers of SNPs (why?)

  # scale the genotypes
  W = apply(EUR2,2,scale)

  # WW'/n
  GRM = W%*%t(W)/length(EUR2p)
  dim(GRM)          # number of individuals by number of individuals

## [1] 503 503

  # plot, first get in small subset into 'long' format then plot
  data = melt(GRM[1:20,1:20])
  ggplot(data, aes(Var1, Var2, fill= value)) + 
    geom_tile() +
    scale_fill_gradient(low = "white", high = "red") +
    coord_fixed()

GRM from GCTA .gz format, with some relatives

Mostly, genomic relationship matrices are made with software such as GCTA. Please see equation 3 in the main GCTA paper for details on how the GRM is calculated in this program. Usually GRMs are stored in (non-human readable) binary format but the following is in g-zipped .grm.gz format which can be easily read into R.

There are two paired files, one called families.grm.gz which contains the lower triangle elements of the GRM (remember the GRM is symmetrical) plus the diagonal elements. The columns in the families.grm.gz file are row index, column index, number of SNPs and \(\pi\) or the relationship value. The file is in .gz format which can be viewed in unix using cat file.gz | head for example.

The corresponding id file (familes.grm.id) contains the family ID and individual ID as columns 1 and 2. Individuals require unique individual IDs but may have the same family ID number.

   # read in the GRM
   GRM = read.table(paste0(dir,"families.grm.gz"))  # lower triangle
   head(GRM) ; tail(GRM)                # 629x629 square matrix, lower triangle provided

##   V1 V2     V3           V4
## 1  1  1 316701  0.999557900
## 2  2  1 316701 -0.004836143
## 3  2  2 316701  1.006090000
## 4  3  1 316701  0.003647990
## 5  3  2 316701 -0.003971321
## 6  3  3 316701  1.006909000

##         V1  V2     V3            V4
## 198130 629 624 316701 -0.0064965470
## 198131 629 625 316701 -0.0039963370
## 198132 629 626 316701 -0.0042707180
## 198133 629 627 316701 -0.0001947815
## 198134 629 628 316701  0.0047688470
## 198135 629 629 316701  0.9967397000

   IDs = read.table(paste0(dir,"families.grm.id"))  # ids, n = 629 matching .grm.gz file
   head(IDs)

##   V1 V2
## 1 18 18
## 2 29 29
## 3 45 45
## 4 50 50
## 5 56 56
## 6 57 57

   dim(IDs)

## [1] 629   2

We are going to take a look at the GRM, and plot the histogram of the diagonal and off-diagonal elements. We can see that this GRM has a few close relatives (larger values in plot below which are not on the diagonal), but most of the individuals are what is usually defined as ‘unrelated’ (\(\pi < 0.05\)) in human genetics.

Note that:
* mean of the diagonal elements, is approx. 1
* mean of the off-diagonal elements, is approx. 0

For a GWAS study, there are two options for handling relatives. Either remove one member of the relative pair (probably appropriate in this case as there are only 20 pairs of relatives) or a mixed model approach could be used. We look at using both of these approaches in the gwas prac session.

  #  small subset & name columns to plot
   names(GRM) = c("ID1","ID2","nSNPs","value")
   GRM2 = GRM[GRM$ID1<101&GRM$ID2<101,]
   ggplot(GRM2, aes(ID1, ID2, fill= value)) + 
      geom_tile() +
      scale_fill_gradient(low = "white", high = "red") +
      coord_fixed()

   # diagonal elements, relationship each individual with themselves
   diags = GRM[GRM$ID1==GRM$ID2,"value"]
   mean(diags)                              # mean diagonals ~ 1

## [1] 1.001292

   hist(diags,breaks=50)

   # off-diagonal elements, relationship each individual with each other
   offdiags = GRM[GRM$ID1!=GRM$ID2,"value"]
   mean(offdiags)                           # mean off-diag ~ 0

## [1] 0.0001750914

   hist(offdiags,breaks=50)

   length(offdiags[offdiags>0.05])          # there are about 400 relative pairs, pi>0.05

## [1] 400

   hist(offdiags[offdiags>0.05],breaks=50)  # looks like most pairs are mostly cousins or less, why?

Part 4: PCA Analysis

PCA analysis is an approach often used to identify population structure in genomic data. You will be provide with a pre-computed genomic relationship matrix for a small subset of the 1000 Genomes project data. Further information on the 1000 Genomes Project can be found here.

The genomic relationship matrix (GRM) in this example has been computed for 50 individuals identified from one of 5 ancestry groups. It was computed from about 1.3M SNP spread evenly throughout the human genome using standard methods in GCTA.

If you need an introduction in Principal Component Analysis (PCA) there is a great 20minute youtube clip called PCA… Step-by-Step by StatQuest. StatQuest also has many other useful bite-size tutorials on topics related to genetics. They are very helpful!

In short, PCA looks for the dimension in your data which explains the most variation in your data. It defines this direction using ‘PC1’ or the first eigenvector. It also returns the amount of variation explained by the eigenvector as an eigenvalue 1. After the first eigenvector/eigenvalue is defined, it looks for the next uncorrelated direction which explains the most variation, called PC2. And so on. Hence a PCA analysis returns pairs of eigenvectors (i.e. the direction of variation) and eigenvalues (related to the variance explained by the eigenvector) which represent the main sources of variation in the original matrix. Formally, it is a decomposition technique used to simplify large matrices. Thus, a square symmetric matrix \(\mathbf{A}\) can be decomposed as:

\[ \mathbf{A} = \mathbf{Q} \boldsymbol{\Lambda} \mathbf{Q}^\top \] where \(\mathbf{Q}\) is the matrix of eigenvectors (columns) and \(\boldsymbol{\Lambda}\) is the diagonal matrix of eigenvalues. The vectors in Q are orthogonal (meaning uncorrelated) with each other while the sum of eigenvalues as a proportion of the sum of eigenvalues gives you the proportion of variation explained by the corresponding vector. It is easy enough to do a PCA in R for a small dataset.

In GWAS, the first few principle components from the GRM are often used as covariates in the analysis to account for population structure. We’ll now look at how to calculate PCs from a pre-computed GRM.

The files used in this part of the prac are:

1kg_winterSchool.grm.gz → pre-made genomic relationship matrix for the 1000G data.

1kg_winterSchool_superPop.txt → file defining the continental populations of the 1000G samples.

First we will read in the data, and have a quick look at the format of the file.

  # read in the data
  x=read.table(paste0(dir,"1kg_winterSchool.grm.gz"))
  labels=read.table(paste0(dir,"1kg_winterSchool_superPop.txt"))[,1]

  # columns are row number, column number, number of SNP, genomic relationship
  head(x)

##   V1 V2      V3         V4
## 1  1  1 1365790 0.95116320
## 2  2  1 1365790 0.06590147
## 3  2  2 1365790 0.92498110
## 4  3  1 1365790 0.05928111
## 5  3  2 1365790 0.07948754
## 6  3  3 1365790 0.92168420

  tail(x)

##      V1 V2      V3          V4
## 1270 50 45 1365790 0.001001333
## 1271 50 46 1365790 0.046047610
## 1272 50 47 1365790 0.034500270
## 1273 50 48 1365790 0.049708950
## 1274 50 49 1365790 0.048485100
## 1275 50 50 1365790 0.878313700

The GRM is a square symmetrical matrix. The data read in are the lower triangular and diagonal elements, as indicated by the row and column number. How many elements do you expect in this file for a 50x50 GRM?

In the first column there will be \(n\) elements, \(n-1\) in the 2nd column, \(n-2\) in the 3rd, etc. until \(1\) element in the 50th column. Check this calculation against the dimensions of the GRM provided.

  # n = number of individuals
  n=50

  # expected number of elements in lower-triangular matrix, inc. diagonals
  N=sum(n:1) ; N

## [1] 1275

  # how many elements read in by R for the lower-triangle + diagonal.
  dim(x)

## [1] 1275    4

Each of the 1000G participants has a ancestry label assigned to it. Use the table() function to count the labels, note that AFR = African, AMR = admixed American, EAS = east Asian, SAS = south Asian, EUR = European.

  table(labels)

## labels
## AFR AMR EAS EUR SAS 
##   7   8  12  11  12

Now construct the GRM as a square matrix so that we can use the eigen() function in R to do the eigenvalue decomposition. The relationship values are in the 4th column of x. We can access the 4th column by using the notation object[rows,columns], where x[,4] will access all rows of the 4th column or x[1,4] will access a single element from the 1st row and 4th column.

  # define the matrix using NA, then fill the matrix
  GRM=matrix(NA,nrow=n,ncol=n)
  for (k in 1:N) {
    i=x[k,1] ; j=x[k,2]
    GRM[i,j]=x[k,4]
    GRM[j,i]=x[k,4]
  }

Now conduct the PCA analysis using the eigen() function. Save the output and check the structure of the returned object using str(). If the matrix is full-rank, we expect 50 eigenvalues with values > 0, and 50 corresponding vectors of length 50 (i.e. a 50x50 matrix). We’ll plot the first 2 PCs with the labels provided and detailing how much variation the PCs capture.

  #eigen value decomposition
   PCA=eigen(GRM)

  #eigenvalues (vector), eigenvectors (50x50 matrix, 1 column per eigenvector)
  str(PCA)

## List of 2
##  $ values : num [1:50] 4.78 2.97 1.46 1.43 1.39 ...
##  $ vectors: num [1:50, 1:50] 0.0406 0.0471 0.049 0.0482 0.1006 ...
##  - attr(*, "class")= chr "eigen"

  # proportion of variation explained by PCs
   total=sum(PCA$values)
   plot(PCA$values[1:10]/total,
        las=1, type="b", lwd=1.2, col="blue", xlab="principal component",
        ylab="Proportion of total variance", pch=20)

  # plot 
  labels=as.factor(labels)
  PC1=PCA$vectors[,1] ; PC2=PCA$vectors[,2]
  plot(PC2~PC1, col=labels, pch=20, las=1, 
       xlab="PC1, 9.2% of variation", ylab="PC2, 5.7% of variation")
  legend("bottomleft",legend=levels(labels),fill=1:5)

In practice, we would use software such as GCTA or PLINK at the command line to (1) calculate the GRM, and (2) conduct the PCA analysis. These programs are much more efficient compared to R. However, we have done it in R here as it is useful for understanding what’s going on.

UQWS Genetic Mapping, session 2 - study design

Kathryn Kemper

30/6/2026