gemma_analysis.Rmd

---
title: "Gemma analysis"
author: "Johan Zicola"
date: "`r Sys.Date()`"
output: html_notebook
---

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

```{r}

# List of packages required for this analysis
pkg <- c("qqman", "ggplot2")

# Check if packages are not installed and assign the
# names of the packages not installed to the variable new.pkg
new.pkg <- pkg[!(pkg %in% installed.packages())]

# If there are any packages in the list that aren't installed,
# install them
if (length(new.pkg)) {
  install.packages(new.pkg, repos = "http://cran.rstudio.com")
}


#Library to plot Manhattan plots
library(qqman)

library(ggplot2)


```

# Function to plot Manhattan plot

The function "GWAS_run.R" modifies the original 'manhattan' function of qqman package by generating
a Bonferroni threshold (5% alpha risk / nb of SNPs analyzed). It will also retrieve
all SNPs above this threshold as output of the function

```{r}
source("GWAS_run.R")

```


###Load output from Gemma
```{r}

dir_file <- "T:/dep_coupland/grp_hancock/johan/GWAS/output/"

file.name <- "Area_nucleus.assoc.clean.txt"

path.file <- paste(dir_file, file.name, sep="")

gwas.results <- read.delim(path.file, sep="\t")


```

###QQ plot of the p-values
```{r echo=FALSE}
qq(gwas.results$P, main=file.name)
```

```{r echo=FALSE}
plot(-log(gwas.results$P)~gwas.results$CHR, main=file.name)
```

```{r echo=FALSE}
hist(-log(gwas.results$P), main=file.name)
```

### Number of SNP per chromosome
```{r echo=FALSE}
as.data.frame(table(gwas.results$CHR))
```

### Manhattan plot

Get a vector of the SNPs with significant value and display the Manhattan plot with bonferroni
 threshold and the SNP "1:24018" highlighted in green

```{r echo=FALSE}
SNP_significant <- GWAS_run(path.file, threshold_pvalue = "bonferroni", highlighted_SNP="1:24018")
```