---
title: "Exercise 2 SISG Association Mapping"
author: "null"
date: "2/6/2017"
output: html_document
---

```{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)
```


# Load the bpdata.csv data file into your R session.  

LHON.txt data file into your R session.  

Can read the data directly from the website if your computer is connected online 

```{r}
bpdata=read.csv("http://faculty.washington.edu/tathornt/sisg/bpdata.csv",header=TRUE)

```

 If file is saved on your computer, could also read the data in from the directory that contains the file.
 
```{r}
bpdata2=read.csv("bpdata.csv",header=TRUE)
```

View the first few lines of the bpdata
```{r}
head(bpdata)
```


Get information about the types of variables in the LHON data frame

```{r}
str(bpdata)
```

#1 Perform a linear regression in R of systolic blood pressure (sbp) on SNP3 using lm() for different models.  


## Additive model

First create a genotype variable with an additive coding based on the counts of the number of T alleles 
```{r}
bpdata$n.minor <- (bpdata$snp3=="TC") +(bpdata$snp3=="CT")+2*(bpdata$snp3=="TT")

head(bpdata$n.minor)
```

Now perform a linear regression analysis with the additive genotype coding variable

```{r}
lm1 = lm(sbp~n.minor, data=bpdata)
summary(lm1)
```
Can obtain confidence intervals of the effect estimates using *confint.default()* function 
```{r}
confint.default(lm1)
```


Note that we obtain the same inference on the association with an additive model when the genotype variable is coded based on the number of copies of the C allele
```{r}
bpdata$n.minorC<- (bpdata$snp3=="TC") + (bpdata$snp3=="CT")+ 2*(bpdata$snp3=="CC")
head(bpdata$n.minorC)
lm1a = lm(sbp~n.minorC, data=bpdata)
summary(lm1a)

```


# Dominant effect model for T allele.  

First create a variable that has a value of 1 for genotypes with at least 1 T allele, and 0 for the genotype that is CC

```{r}
bpdata$domvar <- (bpdata$snp3=="TC") | (bpdata$snp3=="TT")
with(bpdata, table(domvar, snp3))
```

Now perform the linear regression analysis  for  the dominant coded genotype variable and obtain confidenc intervals.

```{r}

lm2 <- lm(sbp~domvar, data=bpdata)
summary(lm2)
confint.default(lm2)
```


# Recessive Model for T allele
```{r}
bpdata$recvar <- (bpdata$snp3=="TT")
with(bpdata, table(recvar, snp3))
lm3<-lm(sbp~recvar, data=bpdata)
summary(lm3)
confint.default(lm3)
```


# 2 Degree of freedom model

```{r}
lm4<-lm(sbp~snp3, data=bpdata)
summary(lm4)
confint.default(lm4)
```

To test the 2 degree of freedom model, can use ANOVA (analysis of variance), which compares the full model to the null model.   First need to remove missing data, so comparing models on the same dataset

```{r}
clean.data <- bpdata[complete.cases(bpdata[,c("snp3")]),]
lm4<-lm(sbp~snp3, data=clean.data)
lm.null  <- lm(sbp~1, data=clean.data)
anova(lm.null, lm4) 
```







# 2. Provide plots illustrating the relationship between sbp and the three genotypes at SNP3. 

Plot of SBP versus minor allele count for snp3 with linear regression line

```{r}
plot(sbp~n.minor, data=bpdata,col="darkgreen")
abline(lm1, col="red", lwd=2)
```

Using *jitter()* adds noise to the points so that we can better see how many observations there are within each genotype group 

```{r}

plot(sbp~jitter(n.minor), data=bpdata,col="darkgreen")
abline(lm1, col="red", lwd=2)
```

Below is a boxplot with a regression line

```{r}

boxplot(sbp~n.minor,data=bpdata, col=c("purple","green","darkblue"))
legend("topright",c("CC","CT","TT"),fill=c("purple","green","darkblue"))
abline(lm1, lwd=2,col="red")

```






## 3. Now redo the linear regression analysis of sbp from question 2 for the additive model, but this time adjust for sex and bmi.   Do the results change?  
```{r}
lm1adj = lm(sbp~n.minor+sex+bmi, data=bpdata)
summary(lm1adj)
confint.default(lm1adj)
```


# 4.  Perform a linear regression for additive model using all SNPs


```{r}
lmall<- lm(sbp~snp1+snp2+snp3+snp4+snp5+snp6+snp7+snp8+snp9+snp10+snp11,data=bpdata)
summary(lmall)
```


Below is the model above with sex and bmi included  as predictors

```{r}
lmalladj<- lm(sbp~snp1+snp2+snp3+snp4+snp5+snp6+snp7+snp8+snp9+snp10+snp11+sex+bmi,data=bpdata)
summary(lmalladj)
```


Below is Another way to do the regression analysis in R without having to write all of the predictors:

```{r}
lm5<- lm(sbp~.,data=bpdata[,c(2:3,5:16)])
summary(lm5)
```




