---
title: "Exercise 6 SISG Association Mapping: Get GRM and Null Model for Transferrin"
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
)
```


First, make sure that you have installed the Bioconductor R packages gdsfmt, SNPRelate, GWASTools, and GENESIS : 
```{r}
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(c("GWASTools", "gdsfmt","SNPRelate","GENESIS"))

```




Load the packages gdsfmt, SNPRelate, GWASTools, and GENESIS into your R session 
```{r}
library(gdsfmt)
library(SNPRelate)
library(GWASTools)
library(GENESIS)
```



## NOW we will use functions in the gdsfmt, SNPRelate packages, and GENESIS to perform a linear mixed model association analysis.   Tutorial for the gdsfmt and SNPRelate package can be found here: <http://corearray.sourceforge.net/tutorials/SNPRelate/>.  Tutorials  for performing mixed model association testing with the GENESIS package is here:
<http://bioconductor.org/packages/release/bioc/vignettes/GENESIS/inst/doc/assoc_test.html>.
 



###  Using the SNPRelate package, create a GDS file for R from the PLINK files
```{r}
bedfile<-"Transferrin.bed" 
famfile<-"Transferrin.fam" 
bimfile<-"Transferrin.bim" 
```


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

Open the GDS file and

```{r}
genofile <- snpgdsOpen('Transferrin.gds')

```

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


Read in the phenotype transferrin phenotype file and identify individuals who do not have missing phenotypes


```{r}
TPHENO<-read.table(file="Tr.pheno", header=FALSE)

# Give names to PHENO variables  
 
names(TPHENO)=c("FAMID", "ID", "Transferrin")

# make a data.frame
 
mydat <- data.frame(sample.id = TPHENO$ID, famID=TPHENO$FAMID,transferrin = TPHENO$Transferrin) 
 
# get a list of IDS for individuals who have non-missing transferrin phenotype values
finaldat<-mydat[!is.na(mydat$transferrin),]
keep<-finaldat$sample.id

```

Calculate the genetic relatedness matrix (GRM) for all individuals who have transferrin measurements. Note, that calculating the GRM will take a long time as there are 2,362 individuals with transferrin measurements.  I have provided the grm in a file, but below are the commands I used.   Uncomment the commands to run them on your computer.

```{r}



### Create the GRM

# grm=snpgdsGRM(genofile, sample.id=keep, snp.id=NULL, autosome.only=TRUE, remove.monosnp=TRUE, maf=.05, missing.rate=.05, method="GCTA", num.thread=1L, with.id=TRUE, verbose=TRUE)


### Save the GRM to a file

# save(grm,file='Transferrin_grm.Rdata')
```

Get the Null Model for Transferrin

```{r}

# Load GRM

load('Transferrin_grm.Rdata')
GRM=grm$grm
row.names(GRM)=grm$sample.id
colnames(GRM)=grm$sample.id

annot <- AnnotatedDataFrame(data=finaldat)

#See the annotated data.frame with pData
head(pData(annot))

## Fit the null model

LMMnullmodNEW=fitNullModel(annot,outcome="transferrin",cov.mat=GRM, family=gaussian) 

## Save the null model

save(LMMnullmodNEW, file='Transferrin_standard_lmm_null_model_NEW.Rdata')

# Close gds file

closefn.gds(genofile)

```




