---
title: "Exercise 5 SISG Association Mapping"
author: ''
date: ''
output:
  html_document: default
  pdf_document: default
---

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




```{r echo=FALSE}
knitr::opts_chunk$set(
  message = FALSE,
  warning = FALSE,
  error = FALSE, 
  collapse = TRUE,
  comment = "#",
  fig.height = 4,
  fig.width = 8,
  fig.align = "center",
  cache = FALSE
)
```

```{r,include=FALSE}
library(dplyr)
library(broom)
library(tidyr)
library(ggplot2)
```



First, make sure that you have installed the Bioconductor R packages gdsfmt and SNPRelate: 

* source("http://bioconductor.org/biocLite.R")
* biocLite("gdsfmt")
* biocLite("SNPRelate")


Load the pakcages gdsfmt and SNPRelate into your R session 
```{r}
library(gdsfmt)
library(SNPRelate)
```




### Investigate the contents of the PLINK files .fam and .bim files

Look at the .fam file and obtain the number of individual in the study
```{r}
FAM<-read.table(file="../SISGAssocData/YRI_CEU_ASW_MEX_NAM.fam",sep=" ", header=FALSE,na="NA")
head(FAM)
dim(FAM)
unique(FAM$V1)
```


Now obtain information about the number SNPs in the .bim file 
```{r}
map<-read.table(file="../SISGAssocData/YRI_CEU_ASW_MEX_NAM.bim",sep="\t", header=FALSE,na="NA")
head(map)
dim(map)
```




### Read in file with Information on Populations of Sample Individuals
 
 
```{r}
POPINFO=read.table(file="../SISGAssocData/Population_Sample_Info.txt",header=TRUE)
```

Obtain a table with the number of sample individuals from each population
```{r}
table(POPINFO$Population)
```


Check to make sure the order of individuals is the same as the PLINK file


```{r}
identical(POPINFO$IID,FAM$V2)
```


## NOW we will use functions in the gdsfmt and SNPRelate packages to perform a PCA.   A tutorial on gdsfmt and SNPRelate can be found at: http://corearray.sourceforge.net/tutorials/SNPRelate/


###  Create a GDS file for R from the PLINK files
```{r}
bedfile<-"../SISGAssocData/YRI_CEU_ASW_MEX_NAM.bed" 
famfile<-"../SISGAssocData/YRI_CEU_ASW_MEX_NAM.fam" 
bimfile<-"../SISGAssocData/YRI_CEU_ASW_MEX_NAM.bim" 
```




### Convert PLINK FILES INTO GDS FORMAT FOR R
```{r}
snpgdsBED2GDS(bedfile, famfile, bimfile, "IPCgeno.gds")
```


Open the GDS file and

```{r}
genofile <- snpgdsOpen("IPCgeno.gds")
```

Explore the GDS file and a few of the components 
```{r}
head(genofile)
head(read.gdsn(index.gdsn(genofile, "sample.id")))
head(read.gdsn(index.gdsn(genofile, "snp.id")))
```

Extract genotype for a specific SNP and plot a histrogram of the genotype values 
```{r}
g <- snpgdsGetGeno(genofile, snp.id="rs3813199")
hist(g)
```



## Now perform a PCA using a function from the SNPRelate package
```{r}
pca <- snpgdsPCA(genofile)
```

## Extract the first 2 PCs and make a data frame

```{r}
tab <- data.frame(sample.id = pca$sample.id,
    EV1 = pca$eigenvect[,1],    # the first eigenvector
    EV2 = pca$eigenvect[,2],    # the second eigenvector
    stringsAsFactors = FALSE)
head(tab)
```

## Make a plot of the top 2 PCs

```{r}
plot(tab$EV2, tab$EV1, xlab="eigenvector 2", ylab="eigenvector 1")
```



### Now incorporate population membership for the PCA plot

Get sample id and create a new dataframe with the population membership
```{r}
sample.id <- read.gdsn(index.gdsn(genofile, "sample.id"))
population=POPINFO$Population
tab <- data.frame(sample.id = pca$sample.id,
    pop = factor(population)[match(pca$sample.id, sample.id)],
    EV1 = pca$eigenvect[,1],    # the first eigenvector
    EV2 = pca$eigenvect[,2],    # the second eigenvector
    stringsAsFactors = FALSE)
head(tab)
```

Make a plot of the top 2 PCs, with points color coded by population membership 
```{r}
plot(tab$EV2, tab$EV1, col=as.integer(tab$pop), xlab="eigenvector 2", ylab="eigenvector 1", main="PCA using all SNPs")
legend("topleft", legend=levels(tab$pop), pch="o", col=1:(nlevels(tab$pop)))
```




### Now make scatterplots of the top 4 PCs with proportional variance explained included

```{r}
pc.percent <- pca$varprop*100
head(round(pc.percent, 2))
lbls <- paste("PC", 1:4, "\n", format(pc.percent[1:4], digits=2), "%", sep="")
pairs(pca$eigenvect[,1:4], col=tab$pop, labels=lbls)


```



### Do we actually need 150K SNPs for population structure infererence in this sample?

Obtain a subset of SNPs based on linkage disequilibrium (LD threshold of 0.2 


```{r}
snpset <- snpgdsLDpruning(genofile, ld.threshold=0.2)
snpset.id <- unlist(snpset)
```


Now perform a PCA with the subset of SNPs
```{r}
pca2 <- snpgdsPCA(genofile, snp.id=snpset.id)
tab2 <- data.frame(sample.id = pca2$sample.id,
    pop = factor(population)[match(pca2$sample.id, sample.id)],
    EV1 = pca2$eigenvect[,1],    # the first eigenvector
    EV2 = pca2$eigenvect[,2],    # the second eigenvector
    stringsAsFactors = FALSE)
```


Make a plot of the top 2 PCs with the subset of SNP, with points color coded by population membership 

```{r}
plot(tab2$EV2, tab2$EV1, col=as.integer(tab$pop), xlab="eigenvector 2", ylab="eigenvector 1", main="PCA using a subset of SNPs")
legend("bottomright", legend=levels(tab$pop), pch="o", col=1:(nlevels(tab$pop)))
```



# Estimate proportional ancestry from PCA

## First Estimate Native American and European Ancestry for the MXL from the PCA.  For simplicity assume that MXL have negligible African Ancestry.


Obtain average PC 2 value for CEU (European) and NAM (Native Americans), since PC2 is relecting European and Native American ancesty
```{r}

avgCEU2=mean(pca2$eigenvect[population=="CEU",2])
avgNAM2=mean(pca2$eigenvect[population=="NAM",2])

```

Estimate proportional MXL ancestry, based on their relative value to the average CEU and NAM PC 2 values

```{r}
MXLadmix=(pca2$eigenvect[population=="MXL",2]-avgNAM2)/(avgCEU2-avgNAM2)
```


Make a table of the estimated admixture proportions 
```{r}
tab2=cbind(MXLadmix,1-MXLadmix)
myorder=order(MXLadmix)
temp=t(as.matrix(tab2[myorder,]))
```


MAKE A BARPLOT OF MXL  ESTIMATED ANCESTRY FROM THE PCA 

```{r}
barplot(temp, col=c("blue","green"),xlab="Individual ", ylab="Ancestry", border=NA,axisnames=FALSE,main="Ancestry of MXL",ylim=c(0,1))
legend("bottomright", c("European","Native American"), lwd=4, col=c("blue","green"), bg="white",cex=0.85)

```

## How about predicting ancestry for 3 ancestral populations? 

See documentation about this here: https://www.rdocumentation.org/packages/SNPRelate/versions/1.6.4/topics/snpgdsEIGMIX


##First define groups
```{r}
groups <- list(CEU = sample.id[population == "CEU"],
               YRI = sample.id[population == "YRI"],
               NAM = sample.id[population=="NAM"])
```

### Can use the snpgdsAdmixProp function to estimate the proportional ancestry from the pca that we did. See
```{r}
prop <- snpgdsAdmixProp(pca2, groups=groups,bound=TRUE)
plot(prop[, "YRI"], prop[, "CEU"], col=as.integer(tab$pop),
     xlab = "Admixture Proportion from YRI",
     ylab = "Admixture Proportion from CEU")
abline(v=0, col="gray25", lty=2)
abline(h=0, col="gray25", lty=2)
abline(a=1, b=-1, col="gray25", lty=2)
legend("topright", legend=levels(tab$pop), pch="o", col=1:5)


```




### Can also obtain barplots of proportional ancestry of  MXL  that was estimated from PCA ###

```{r}
MXLprop<-data.frame(prop[population=="MXL",])
myorder=order(MXLprop$CEU)
temp=t(as.matrix(MXLprop[myorder,]))


barplot(temp, col=c("blue","red","green"),xlab="Individual ", ylab="Ancestry", border=NA,axisnames=FALSE,
        main="Estimated Ancestry of MXL from PCA",ylim=c(0,1))
legend("bottomright", c("European","African","Native American"),lwd=4, col=c("blue","red","green"),bg="white",cex=0.85)

```






