---
title: "Exercise 6 SISG Association Mapping: Run LMM 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
)
```


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

```{r}

load(file='Transferrin_standard_lmm_null_model_NEW.Rdata')

# Read in gds file 
geno <- GdsGenotypeReader(filename = 'Transferrin.gds')

##
TPHENO<-read.table(file="Tr.pheno", header=FALSE)

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


# make a scan annotation data.frame
mydat <- data.frame(scanID = TPHENO$ID, famID=TPHENO$FAMID,transferrin = TPHENO$Transferrin)
scanAnnotTRANS=ScanAnnotationDataFrame(mydat)



# create a GenotypeData class object with paired ScanAnnotationDataFrame

genoData <- GenotypeData(geno, scanAnnot = scanAnnotTRANS)



# Create a genotype iterator for testing chromosome 22


chr <- getChromosome(genoData)
snpID <- getSnpID(genoData)
chr22 <- snpID[chr == 22]
genoIterator <- GenotypeBlockIterator(genoData, snpInclude=chr22)

#Start Association test for chromosome 22 SNPs
assoc <- assocTestSingle(genoIterator, null.model = LMMnullmodNEW, test = "Wald")

TOP <- assoc[order(assoc$Wald.pval),] 
head(TOP,10,)


#### No do associaton analysis fo all of the genome.  This takes a long time  depending on your computer


genoIterator <- GenotypeBlockIterator(genoData)

#########################################################

#Start Association test for all SNPs
assoc <- assocTestSingle(genoIterator, null.model = LMMnullmodNEW, test = "Wald")


write.csv(assoc,file="LMM_Transferrin_Assoc_Results_NEW.csv",row.names=F,quote=F)


TOP <- assoc[order(assoc$Wald.pval),] 
head(TOP,10,)

# close the genotype data file
close(genodata)


```





