---
title: 'Principal Component Analysis of Genotypic Data'
author: "International Statical Genetics Workshop"
#date: '2023-03-07'
output:
  html_document: default
  pdf_document: default
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

<span style="color:red;">**This practical runs on the workshop RStudio server.**</span>

<!-- This practical aims to strengthen your understanding of the effect of genetic drift **Lecture 6** and how it can lead to population stratification.  -->

First login to RStudio

```{bash, eval=FALSE}
https://workshop.colorado.edu/rstudio/
```

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

```{R, eval=TRUE}
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`

```{R, eval=TRUE}
system("ls -lt")
```

or directly via the terminal
```{bash, eval=TRUE}
ls -lt
```


We now read the labels for each individual in the data.
```{R, eval=TRUE}
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. 

```{R, eval=TRUE}
system("plink --bfile 1kg_hm3_qced --pca 10 --out myFirstPCA")
```

Now let's read and plot the results in `R`
```{R, eval=TRUE}
## 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)
```

<span style="color:blue;">
*Question 1a. How many (ancestry) clusters do you see?*
</span>

<span style="color:blue;">
*Question 1b. What is the ancestry of someone with a PC1 value of 0? Same question for a PC2 value 0.*
</span>

Let's use k-means approach to determine the clusters
```{R, eval=TRUE}
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...
```{R, eval=TRUE}
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)
```


<span style="color:blue;">
*Question 2. Are the k-means clusters concordant with geographical labels? If not then how can you improve concordance?*
</span>

To answer the question you can run the `R` code below

```{R, eval=TRUE}
## 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...
```{R, eval=TRUE}
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)
  }
}
```

So let's redefine k-mean clusters with more PCs...
```{R, eval=TRUE}
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"])
```

<span style="color:blue;">
*Question 3. Can we improve concordance between PC-based clusters and geographical labels when using more than 5 PCs?*
</span>

<span style="color:blue;">
*Question 4. What can you conclude overall?*
</span>

<!-- 1) PCA captures geographical information (mainly between continents) -->
<!-- 2) The more dimensions the better but don't expect 100% with geographical labels. -->

## 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
```{R, eval=TRUE}
## 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...

```{R, eval=TRUE}
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...
```{R, eval=TRUE}
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`

```{R, eval=TRUE}
system("head myPCA_noAMR_snp_loading.pcl")
```
or the terminal
```{bash, eval=TRUE}
head myPCA_noAMR_snp_loading.pcl
```

Step 3. Project AMR samples...
```{R, eval=TRUE}
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`

```{R, eval=TRUE}
system("head AMR_PCA.proj.eigenvec")
```

or using the terminal

```{bash, eval=TRUE}
head AMR_PCA.proj.eigenvec
```

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

```{R, eval=TRUE}
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)

```

<span style="color:blue;">
*Question 5. Compare PC coordinates of AMR samples when they are included versus excluded from PCA*
</span>

```{R, eval=TRUE}
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"])

```


<span style="color:blue;">
*Question 6. Assign an ancestry cluster to the AMR samples*
</span>

```{R, eval=TRUE}
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)
```