Principal Component Analysis of Genotypic Data

This practical runs on the workshop RStudio server.

First login to RStudio

https://workshop.colorado.edu/rstudio/

Copy the Practical document into your home directory using the following R commands.

setwd("~/.")
system("mkdir -p Day1")
system("cp /home/practicals/1.4.PrincipalComponentAnalysis_LoicYengo/final/Practical_PCA.html ~/Day1/.")

## Now copy the data
system("mkdir -p Day1/pca_prac")
system("cp /home/practicals/1.4.PrincipalComponentAnalysis_LoicYengo/final/data/1kg_* ~/Day1/pca_prac/.")
system("cd Day1/pca_prac")
setwd("Day1/pca_prac")

Data description

In this practical we use data from the 1000 Genomes Project (Phase 3). The dataset consists of 2,504 individuals from 26 diverse populations grouped into 5 super-population (AFR: African; AMR: Admixed American; EAS: East-Asian; EUR: European; SAS: South-Asian). More details about the 1000 Genomes data can be found at https://www.internationalgenome.org/. For the practical, we will QC genotypes for 80,244 SNPs.

Check that all the files are copied correctly using R

system("ls -lt")

We now read the labels for each individual in the data.

labels <- read.table("1kg_pop2504.txt",h=T)
head(labels)

Part 1: PCA with PLINK and ancestry calling

Run the following PLINK command (called from R, i.e., ghe actual command is everything between system(\" PLINK COMMAND \")) to calculate 10 PCs from the set of genotypes stored in the 1kg_hm3_qced.*files.

system("plink --bfile 1kg_hm3_qced --pca 10 --out myFirstPCA")

Now let’s read and plot the results in R

## Read the results
pca_res <- read.table("myFirstPCA.eigenvec")
colnames(pca_res) <- c("FID","IID",paste0("PC",1:10))

## Read eigen values
eigenvals <- read.table("myFirstPCA.eigenval")[,1]

## Check the results
head(pca_res)
npc <- 10
for(k in 1:npc){
  pca_res[,paste0("PC",k)] <- pca_res[,paste0("PC",k)] * sqrt(eigenvals[k])
}

## Add labels
pca_res <- merge(pca_res,labels,by.x="IID",by.y="sample")

plot(pca_res[,"PC1"],pca_res[,"PC2"],pch=19,axes=FALSE,xlab="PC1",ylab="PC2",col=1,main="")
axis(1);axis(2)
abline(h=0,v=0,col="grey",lty=2)

Question 1a. How many (ancestry) clusters do you see?

Question 1b. What is the ancestry of someone with a PC1 value of 0? Same question for a PC2 value 0.

Let’s use k-means approach to determine the clusters

ncluster   <- 5
cluster1_2 <- kmeans(pca_res[,c("PC1","PC2")],ncluster,nstart = 10)$cluster
plot(pca_res[,"PC1"],pca_res[,"PC2"],pch=19,axes=FALSE,xlab="PC1",ylab="PC2",col=cluster1_2,main="")
axis(1);axis(2)
abline(h=0,v=0,col="grey",lty=2)

Now let’s use the geographic labels…

Colors  <- c("dodgerblue","coral1","goldenrod","seagreen4","purple3")
Pops    <- names(table(labels$super_pop))
names(Colors) <- Pops

plot(pca_res[,"PC1"],pca_res[,"PC2"],pch=19,axes=FALSE,xlab="PC1",ylab="PC2",col=Colors[pca_res[,"super_pop"]],main="")
axis(1);axis(2)
abline(h=0,v=0,col="grey",lty=2)
legend("topright",legend=Pops,fill=Colors,box.lty=0)

Question 2. Are the k-means clusters concordant with geographical labels? If not then how can you improve concordance?

To answer the question you can run the R code below

## Credit to Museful (stackoverflow)
permutations <- function(n){
  if(n==1){
    return(matrix(1))
  } else {
    sp <- permutations(n-1)
    p <- nrow(sp)
    A <- matrix(nrow=n*p,ncol=n)
    for(i in 1:n){
      A[(i-1)*p+1:p,] <- cbind(i,sp+(sp>=i))
    }
    return(A)
  }
}

concordance <- function(x,y){
  tab <- table(x,y)
  n   <- nrow(tab)
  m   <- ncol(tab)
  if(n != m){
    cat("Error: Number of clusters must be the same for x and y.")
    return(NULL)
  }
  perms <- permutations(n)
  conc  <- sapply(1:nrow(perms),function(i) sum( diag(tab[perms[i,],]) ) ) / sum(tab)
  return(max(conc))
}

concordance(x=cluster1_2,y=pca_res[,"super_pop"])

Let’s now visualize for PCs…

op <- par(mfrow=c(5,2))
for(j in 1:4){
  for(k in (j+1):5){
    plot(pca_res[,paste0("PC",j)],pca_res[,paste0("PC",k)],pch=19,axes=FALSE,
         xlab=paste0("PC",j),ylab=paste0("PC",k),col=Colors[pca_res[,"super_pop"]],
         main=paste0("PC",j," vs PC",k))
    axis(1);axis(2)
    abline(h=0,v=0,col="grey",lty=2)
  }
}
par(op)

So let’s redefine k-mean clusters with more PCs…

ncluster    <- 5
cluster1_5 <- kmeans(pca_res[,paste0("PC",1:5)],ncluster,nstart = 10)$cluster
plot(pca_res[,"PC1"],pca_res[,"PC2"],pch=19,axes=FALSE,xlab="PC1",ylab="PC2",col=cluster1_5,main="")
axis(1);axis(2)
abline(h=0,v=0,col="grey",lty=2)

## Test concordance again
concordance(x=cluster1_5,y=pca_res[,"super_pop"])

Question 3. Can we improve concordance between PC-based clusters and geographical labels when using more than 5 PCs?

Question 4. What can you conclude overall?

Part 2: Projected PCA using GCTA

In this part of the practical, we will calculate PCs without AMR samples, then project back to the PCA map. We will use the software tool GCTA. This procedure has multiple steps. First, step run PCA in without AMR samples.

We first create a list of ID for AMR samples

## Create a list AMR sample
write.table(pca_res[which(pca_res[,"super_pop"]=="AMR"),c("FID","IID")],"AMR.id",quote=F,col.names=F,row.names=F)

## Check the file
system("head AMR.id")

Step 1. Run PCA without AMR samples…

system("gcta64 --bfile 1kg_hm3_qced --remove AMR.id --make-grm --out myGRM_noAMR")
system("gcta64 --grm myGRM_noAMR --pca 2 --out myPCA_noAMR")

Step 2. Calculate SNP loadings…

system("gcta64 --bfile 1kg_hm3_qced --remove AMR.id --pc-loading myPCA_noAMR --out myPCA_noAMR_snp_loading")

Check that the loading file named myPCA_noAMR_snp_loading.pcl has been created. You can use R

system("head myPCA_noAMR_snp_loading.pcl")

Step 3. Project AMR samples…

system("gcta64 --bfile 1kg_hm3_qced --keep AMR.id --project-loading myPCA_noAMR_snp_loading 2 --out AMR_PCA")

Check that the loading file named AMR_PCA.proj.eigenvec has been created. You can use R

system("head AMR_PCA.proj.eigenvec")

or using the terminal

head AMR_PCA.proj.eigenvec

Let’s now load the results in R and visualize them.

ref <- read.table("myPCA_noAMR.eigenvec",stringsAsFactors=FALSE)
amr <- read.table("AMR_PCA.proj.eigenvec",stringsAsFactors=FALSE)
colnames(ref) <- colnames(amr) <- c("FID","IID","PC1","PC2")

## Read eigenvalues
eigenvals <- read.table("myPCA_noAMR.eigenval")[,1]
npc <- 2
for(k in 1:npc){
  ref[,paste0("PC",k)] <- ref[,paste0("PC",k)] * sqrt(eigenvals[k])
  amr[,paste0("PC",k)] <- amr[,paste0("PC",k)] * sqrt(eigenvals[k])
}

ref <- merge(ref,labels,by.x="IID",by.y="sample")

plot(ref[,"PC1"],ref[,"PC2"],pch=19,axes=FALSE,xlab="PC1",ylab="PC2",col=Colors[ref[,"super_pop"]],main="")
axis(1);axis(2)
abline(h=0,v=0,col="grey",lty=2)
legend("topright",legend=Pops,fill=Colors,box.lty=0)

## Plot AMR projections
points(amr[,"PC1"],amr[,"PC2"],col=Colors["AMR"],pch=17)

Question 5. Compare PC coordinates of AMR samples when they are included versus excluded from PCA

comparison <- merge(pca_res[,c("IID","PC1","PC2")],amr[,c("IID","PC1","PC2")],by="IID")

apply(comparison[,-1],2,mean)
apply(comparison[,-1],2,sd)

cor(comparison[,"PC1.x"],comparison[,"PC1.y"])
cor(comparison[,"PC2.x"],comparison[,"PC2.y"])

Question 6. Assign an ancestry cluster to the AMR samples

clus_centers  <- aggregate(cbind(PC1,PC2)~super_pop,data=ref,FUN=mean)

## Distance calculate distance to each cluster
Distances <- matrix(NA,nrow=nrow(amr),ncol=4)
for(i in 1:nrow(amr)){
  for(k in 1:4){
    Distances[i,k] <- sum( ( as.numeric(amr[i,c("PC1","PC2")]) - as.numeric(clus_centers[k,-1]) )^2 )
  }
}
clus_amr <- clus_centers[apply(Distances,1,which.min),1]
table(clus_amr)

Principal Component Analysis of Genotypic Data

International Statical Genetics Workshop

Data description

Part 1: PCA with PLINK and ancestry calling

Part 2: Projected PCA using GCTA