---
title: "vGWAS basic tutorial"
author:
 - "Xia Shen"
 - "Kristian Ullrich^[Max Planck Institute For Evolutionary Biology, ullrich@evolbio.mpg.de]"
date: "`r Sys.Date()`"
output: rmarkdown::html_vignette
bibliography: vGWAS.bib
vignette: >
  %\VignetteIndexEntry{vGWAS basic tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)
```

## Vignette - vGWAS basic tutorial

### A Brief Tutorial of the R Package vGWAS

## Setup Notes

To use the vGWAS package, of course, an R environment is required. Visit:

http://www.r-project.org

and install R for the operating system.

Start R and in the R console, type the following command to install the package dependencies first:

```{r}
#install.packages("devtools")
#install.packages("knitr")
#install.packages("dglm")
#install.packages("doParallel")
#install.packages("foreach")
#install.packages("genio")
#install.packages("hglm")
#install.packages("Matrix")
#install.packages("onewaytests")
#library(devtools)
#devtools::install_github("kullrich/vGWAS", build_vignettes = TRUE, dependencies = FALSE)
```

Now the package is installed in the R library.

## Example

Load the library in the R console, and open this vignette with the command:
```{r}
library(vGWAS)
```

Nine main functions in the package are ready to use:

1. bfmedian.test
2. brown.forsythe.tes
3. getMAF
4. vGWAS
5. vGWASparallel
6. vGWAS.gc
7. vGWAS.variance
8. plot (S3 method for vGWAS object)
9. summary (S3 method for vGWAS object)

Run the following commands to load the example data:
```{r}
data(pheno)
data(geno)
data(chr)
data(map)
```
pheno is a numeric vector of the simulated phenotypic values. By running:

```{r}
hist(pheno, breaks = 30, density = 15, col = "darkred")
```

The command:
```{r}
table(chr)
```
shows exactly the number of markers on each of the five simulated chromosomes.

To pre-filter the `geno.matrix` data for minor-allele-frequency, use the
`getMAF` function as follows:
```{r}
data(geno)
dim(geno)
maf <- getMAF(geno.matrix = geno, geno.snp = "col")
geno.maf <- geno[, maf > 0.05]
dim(geno.maf)
chr.maf <- chr[maf > 0.05]
map.maf <- map[maf > 0.05]
```

Now, the objects loaded in R are ready for a vGWA scan, which can be done using the single command:
```{r}
vgwa <- vGWAS(
  phenotype = pheno,
  geno.matrix = geno,
  marker.map = map,
  chr.index = chr,
  pB = FALSE)
```
A progress bar will indicate the progress of the scan if pB is set to TRUE.

The same can be run using multiple cores.

__Note:__ using the `vGWASparallel` function one need to specify, if the SNPs
in the `geno.matrix` are at columns (`geno.snp = "col"`)
or at rows (`geno.snp = "row"`)
```{r}
#vgwa <- vGWASparallel(
#  phenotype = pheno,
#  geno.matrix = geno,
#  geno.snp = "col",
#  marker.map = map,
#  chr.index = chr,
#  pB = FALSE,
#  ncores = 2)
```

When the scan is finished, all the output statistics will be returned as a list into the object vgwa, which belongs to the class ‘vGWAS’. Any object that has a structure belonging to class ‘vGWAS’ can be directly passed into S3 method function plot. For instance, simply run the following command, we can plot the results in vgwa:
```{r}
plot(vgwa)
```

which produces Figure 2. There is a clear peak above the Bonferroni corrected threshold (dashed orange line).

Regarding the marker that gave the highest score, the heritability explained by the mean and variance can be split and calculated via:
```{r}
# get marker with lowest p.value
marker.lowest <- vgwa[["p.value"]] == min(vgwa[["p.value"]])
vGWAS.variance(
  phenotype = pheno,
  marker.genotype = geno[, marker.lowest])
```

The output can also be stored if assigning the function call to an object.

## Remarks

The package source and further development information are on the R-Forge project page:
https://r-forge.r-project.org/projects/vgwas/

A reimplementation of `vGWAS` introducing the parallel function `vGWASparallel`
as a devtool package can be found here:
https://github.com/kullrich/vGWAS