---
title: "snpReady: a helper tool to run genomic analysis"
author: Italo Granato and Roberto Fritsche-Neto
date: 05/2017
output: 
   rmarkdown::html_document:
      theme: default
      highlight: pygments
      toc: true
      depth: 3
      css: style.css

vignette: >
  %\VignetteIndexEntry{snpReady-vignette}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

<style>
 body {text-align: justify}
 p {line-height: 2em;}
 </style>

# **Introduction**
  
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;In recent years, advances in genomic sequencing through next generation sequencing have enabled the development of millions of new markers, which have been consistently used in studies of important agronomic traits [(Edwards and Batley, 2010)](#R1.1). Hence, breeders have focused on studies that allow the association of markers with the phenotypic of interest, make predictions of performance or in research involving population studies and diversity analysis.  
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;After generating a large amount of genomic data,  the genomic data must pass through a quality control and imputation of the missing genomic data. Some GS models that take advantage of dimensionality reduction like GBLUP and RKHS needs to construct relationship matrices. Moreover, the understanding of population genetics parameters is important as well.  
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Therefore, there is the need of prepare genomic datasets, in such a way that it can be easily applied in a great range of studies. Hence, we propose proposed snpReady package based on needs of setting datasets ready to run genomic studies in leading genomic applications. Thus, we include in this package the three primary critical needs faced before running genomic analyses: [preparation and quality control of datasets](#QC), [estimation of relationship matrices](#GRM) and [estimation of basic population genetics parameters](#PGS).
  
#  **Quality control** {#QC}
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The function was created with the purpose of recoding and reshape the matrices obtained from different SNP genotyping platforms and let it ready to be used in genomic analyses. Thus, it reshapes, recodes, makes quality control and imputation of missing data in the dataset. It also cleans the map based on the same threshold used in the raw data set.  
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Marker matrix used as input can be organized in two ways. In the *long* format, for each subject, there is an observation of each SNP and its alleles. Thus, if there are $n$ individuals and $p$ markers, the matrix is the order of is $(n \times p) \times 4$ where columns represent samples and SNP identification and one for each allele, in this particular order. In order to illustrate the process, we use a maize data set with 64 lines and 539 SNPs.  
 
```{r}
library(snpReady)
geno <- read.table("http://italo-granato.github.io/geno.txt", header = TRUE, na.strings = "NA")
head(geno)
dim(geno)

```
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Another format that can be used as input is *wide*, where samples are in the row and SNPs in columns. In this case, the matrix is the order of $n \times p$.
```{r, echo = FALSE}
tmp <- data.frame(sample = geno$sample,
                  marker = geno$marker,
                  alleles = paste(geno$allele.1, geno$allele.2, sep = ""), stringsAsFactors = F)
library(reshape2)
x <- acast(data = tmp, sample ~ marker, value.var = "alleles")
x[x == "NANA"] <- NA
```
```{r}
x[1:10,1:5]
```

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The raw data can be coded as the nitrogenous base (A, C, G, and T) or the standard A and B. However, if the data was already recoded this can be set by `base` argument and only quality control is made. Thus, if the base is `FALSE` dataset must be coded as 0, 1 and 2. Missing data should be set as NA.  
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Quality control (QC) for genomic data is based on removing individuals and markers with poor information. In general, for individuals, it can be associated with the amount of missing markers. Hence, samples which do not meet some threshold of missing data can be removed through `sweep.sample`. For markers, QC is based on allele frequency and the amount of missing data. Markers with a low frequency of one of its alleles usually are non-informative and in some situations are considered monomorphic. Therefore, they can be removed. The same trend is applied for missing data. Thus, the QC process is made through `MAF` and `call.rate`.   
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Along with removing non-informative markers, imputation of missing data becomes necessary. We present two types of map-independent imputation. One is based on the Wright equilibrium. For a missing position, we assume that the probability of taking up value is dependent on both allele frequency of the SNP and the level of homozygosity of an individual. Thus:
$$P(x_{ij})=\left\{
                \begin{array}{ll}
                  P(x = 0) = (1 - p_j)^2+ p_j (1 - p_j ) F_i\\
                  P(x = 1) = 2p_j (1 - p_j )-2p_j (1 - p_j)\\
                  P(x=2)= p_j^2+p_j (1 - p_j ) F_i
                \end{array}
              \right.$$
Where $p_i$ is the frequency of the major allele for an SNP $i$, and $F_j$ is the level of homozygosity of an individual $j$ estimate as a proportion of the amount of homozygous loci relative to the total of loci. Another method of imputation currently implemented is based on the mean of each SNP. Each missing position of a SNP $j$ is replaced by its mean. Thus:
$$\bar{p}= 2p_j$$
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;In this section we are going to present a basic usage for `raw.data` in a maize data set. It is composed of 64 inbred lines genotyped with 539 SNPs. First, let's run a basic quality control on this data set.  
  
```{r}
geno.ready <- raw.data(data = as.matrix(geno), frame = "long", base = TRUE, sweep.sample = 0.5, call.rate = 0.95, maf = 0.10, imput = FALSE)
M <- geno.ready$M.clean
M[1:10,1:5]
```
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Above, we recode and clean the dataset according to some quality control parameters. Using a call rate of 0.95 means that it can accept markers with only a maximum of 5% of missing data. MAF of 0.1 means that markers with a minor frequency allele less than 0.1 were removed. For individuals,the threshold used for sweep.sample was 0.5, which means that samples with more than 50% of missing data were removed. 

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Besides the cleaned matrix, a report is also outputted as a summary on how many and what markers were removed by the steps in the quality control.
```{r}
geno.ready$report
```
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;For instance, 67 and 151 markers were removed by MAF and call rate, respectively. Moreover, two individuals were removed. Furthermore, the report shows the SNP  and sample IDs removed by each procedure. It is important highlight that some SNPs may fail in both controls applied. Thus, their identification will appear in both sections. Here, 204 markers were removed, which 14 did not attend either QC applied.

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Now, the quality control can be combined with imputation. When using the imputation, it is needed to choose between the supported methods. 
```{r}
geno.ready2 <- raw.data(data = as.matrix(geno), frame = "long", base = TRUE, sweep.sample = 0.5, call.rate = 0.95, maf = 0.10, imput = TRUE, imput.type = "wright", outfile = "012")
Mwrth <- geno.ready2$M.clean
Mwrth[1:10,1:5]
```
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Regarding the Wright's method, its accuracy is about 62%. Imputation based on mean also have intermediate accuracy. These methods are less accurate than those map-dependent methods, such as the one developed in BEAGLE [(Browning and Browning, 2016)](#R1). However, as described by [Rutkoski et al. (2013)](#R7) with lower missing data rates allowed a mere imputation is enough.

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Concerning the output formats, the `outfile` argument is used to export the cleaned matrix in the suitable format. Default is the count of reference allele such that $AA = 2$, $Aa = 1$, $aa = 0$. Another format is coded as -1, 0, 1, case considering $p_j$ as 0.5. Finally, a special format is suitable for use in the `STRUCTURE` software [(Pritchard et al., 2000)](#R6). To generate this output, is necessary the raw data be coded as nitrogenous bases.
```{r}
geno.readySTR <- raw.data(data = as.matrix(geno), frame = "long", base = TRUE, sweep.sample = 0.5, call.rate = 0.95, maf = 0.10, imput = FALSE, outfile = "structure")
Mstr <- geno.readySTR$M.clean
Mstr[1:10,1:5]
```
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;In this output, each sample is split into two rows, one for each allele. Nitrogenous bases are then recoded to a specific number such that A is 1, C is 2, G is 3 and T is 4 and missing data are assigned as -9. Also, given that `STRUCTURE` can handle missing data, arguments related to imputation are ignored when this output is selected.
  
# **Genomic relationship matrix** {#GRM}
  
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Genomic relationship matrices (GRM) are created through the `G.matrix` function. Genomic prediction models use these matrices, especially in the G-BLUP. Different kinship parametrizations were proposed with the aim of increasing the accuracy of prediction of genomic selection. `G.matrix` function estimates four types of additive and one dominant genomic relationship matrix.   
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The matrix used in the data entry should be coded as allele count, where $AA = 2$, $Aa = 1$, $aa = 0$. However, it accepts continuous values for some SNPs. The second argument is the method to be used to generate the GRM. There are four methods to generate it, three forms of GRM and the Gaussian kernel (GK). Methods currently implemented are the one proposed by [VanRaden (2008)](#R8), two methods proposed by [Yang et al. (2010)](#R9), the UAR (Unified Additive Relationship) and adjusted UAR, and the Gaussian kernel [(Pérez-Elizalde et al., 2015)](#R5).
  
```{r}
G <- G.matrix(M = Mwrth, method = "VanRaden", format = "wide") 
Ga <- G$Ga
Ga[1:5,1:5]

Gd <- G$Gd
Gd[1:5,1:5]

```
Just for the vanRaden method, as shown above, two matrices additive and dominance are generated. Otherwise, only the additive matrix is outputted.
```{r}
G <- G.matrix(M = Mwrth, method = "UAR", format = "wide") 
G[1:5,1:5]
```
Two forms of output are generated, one as a matrix of order $n \times n$.
```{r}
dim(G)
```
Another output is the long format, such that the inverse is used to create a table in a suitable format for `ASREML-R`, where three columns are representing the row, columns and respective value of the lower diagonal matrix. 
```{r}
G <- G.matrix(M = Mwrth, method = "UAR", format = "long") 
head(G)
```

These two forms are suitable to use in many aplications like BGLR [(Pérez and De Los Campos, 2014)](#R4), rrBLUP [(Endelman, 2011)](#R3) and ASREML-R [(Butler et al., 2009)](#R2). 
  
# **Population genetics summary** {#PGS}

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The `popgen` function aims to produce a summary about population genetics parameters. Thus, for any marker locus $j$ with alleles $A_1$ and $A_2$, the function estimates:

* Allele frequencies

$\begin{aligned} 
f(A_1) = p_j = \frac{nA_1}{2N} \\
\\
f(A_2) = q_j = 1 - p_j
\end{aligned}$
<br>

* Minor allele frequency
  
$maf = min(p_j, q_j)$  

* Observed heterozygosity
  
$\begin{aligned} 
H_o=\frac{nH_j}{N}
\end{aligned}$
<br>

* Expected heterozygosity 

$H_e = 2p_jq_j$
<br>

* Nei’s genetic diversity

$DG = 1 - p_j^2 - q_j^2$
<br>

* Polymorphic information content

$PIC = 1-(p_j^2 + q_j^2) - (2p_j^2q_j^2)$
<br>

* Missing rate

$\begin{aligned} 
H_o=\frac{nNA_j}{N}
\end{aligned}$
<br>

* Hardy-Weinberg equilibrium statistic

$\begin{aligned} 
\chi^2=\frac{1}{d}\sum_{k=1}^{3} \frac{(O_k - E_k)^2}{E_k}
\end{aligned}$
<br>

where $nA_1$ is the number of copies of $A_1$ allele in the population, $N$ is the number of individuals, $nH_j$ is the number of heterozygous genotypes (of type $A_1A_2$ or $A_2A_1$) in the locus $j$, $O_k$ is the observed values for the genotypes $0$, $1$ and $2$, $E_k$ is the expected values for $0 = N * (1 - p_j)^2$, $1 = N * 2p_j(1-p_j)$ and $2=N * p_j^2$.
  
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Moreover, for any individual $i$, the function provides estimates of:

* Observed heterozygosity

$\begin{aligned}
H_{o_i} = \frac{nH_i}{m}
\end{aligned}$

* Inbreeding coefficient

$\begin{aligned}
F_i=\frac{O(H_i)-E(H)}{m-E(H)}
\end{aligned}$

where $nH_i$ is the number of heterozigous genotypes (of type $A_1A_2$ or $A_2A_1$) in the individual $i$, $m$ is the number of markers, $O(H_i)$ is the observed homozygosity for individual $i$, $E(H) = \sum_{j} 1-2p_j(1-p_j)$ is the expected homozygosity across all snps.

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; Then, based on the estimates described above,  populational parameters, means and bounders of 
DG, PIC, MAF, Ho, F and S are provided. Besides, measures of variability are estimated:

* Effective populational size

$
\begin{aligned}
Ne =\left(\frac{1}{2\bar{F_i}} \right) N
\end{aligned}
$

* Additive variance due to the allele frequencies

$\begin{aligned}
Va = \sum_{j=1}^{m} 2p_jq_j
\end{aligned}$

* Dominance variance due to the allele frequencies

$\begin{aligned}
Vd = \sum_{j=1}^{m}  (2 p_{j} q_{j})^2
\end{aligned}$

where $\bar{F_i} = \frac{\sum_{i=1}^{N} F_i}{N}$ is the mean of $F_i$

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;In general, data set needs to be coded as allele count and missing data can be accepted. Estimates for the whole population is presented in a list. Thus, there are estimates for `Genotype`, `Markers`, `Population` and `Variability`. 
```{r}
pop.gen <- popgen(M = Mwrth)
head(pop.gen$whole$Markers)
head(pop.gen$whole$Genotypes)
head(pop.gen$whole$Population)
head(pop.gen$whole$Variability)
```

&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;The `popgen` function also allows the assignment of subpopulations to individuals. It estimates the same population genetic parameters for each subpopulation, such as effective population size, components of additive and dominance variance, and endogamy. In our example, let's split the whole population into two subpopulations according to the nitrogen use efficiency (NE), the high and low one. 
```{r}
subgroups <- as.matrix(c(rep("HNE", 10), rep("LNE", 52)))
pop.gen <- popgen(M = Mwrth, subgroups = subgroups)
```
  
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;It produces for each group the same estimates described above. For example, in the HNE group:
  
```{r}
head(pop.gen$bygroup$HNE$Markers)
head(pop.gen$bygroup$HNE$Genotypes)
```
and for LNE group:
```{r}
head(pop.gen$bygroup$LNE$Markers)
head(pop.gen$bygroup$LNE$Genotypes)
```
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Moreover, the separation by groups allows the identification of exclusive/absent and fixed alleles in the assigned sub-populations. Which are in the HNE:
```{r}
head(pop.gen$bygroup$HNE$exclusive)
head(pop.gen$bygroup$HNE$fixed)
```
and in the LNE:
```{r}
head(pop.gen$bygroup$LNE$exclusive)
head(pop.gen$bygroup$LNE$fixed)
```
&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;As can be noted, in the HNE group there are many fixed alleles and no particular alleles. On the other hand, in LNE, the opposite was observed. Hence, for two populations, considering that quality control was made and monomorphic markers were removed, if one allele is fixed in one population, the other will be exclusive in another subpopulation.

In the presence of subgroups, `popgen` also produces the Wright's F statistics in order to measure genetic divergence between populations. $F_{IS}$ measures the deficiency of average heterozygotes in each population, $F_{ST}$ is related to the degree of gene differentiation among populations and $F_{IT}$ is the deficiency of average heterozygotes in subpopulations. The parameters are estimated as follows:

$\begin{aligned}
F_{IS} = 1 - \dfrac{H_I}{H_S}
\end{aligned}$
  
  
$\begin{aligned}
F_{ST} = 1 - \dfrac{H_S}{H_T}
\end{aligned}$
  
  
$\begin{aligned}
F_{IT} = 1 - \dfrac{H_I}{H_T}
\end{aligned}$

where $H_I$ is the weighted average of observed heterozygosity in the subpopulations, $H_S$ is the average expected heterozygosity estimated from each subpopulation and $H_T$ is the expected heterozygosity in the total population. Besides the estimates of F parameters for each marker locus, estimates considering means of $H_I$, $H_T$ and $H_S$ for each population is produced. Thus, There is a comparison considering all population and in a pairwise format which each pair of populations are compared against each other.
  
  
# **References**

#### [1] <span> Edwards, D. and J. Batley. 2010. Plant genome sequencing: applications for crop improvement. Plant Biotechnology Journal 8: 2-9.[Available Here](http://www.sciencedirect.com/science/article/pii/S0958166913006435)</span> {#R1.1}

#### [2] <span>Browning, B.L., and S.R. Browning. 2016. Genotype Imputation with Millions of Reference Samples. **American Journal of Human Genetics** 98(1): 116–126. [Available Here](http://www.cell.com/ajhg/abstract/S0002-9297%2815%2900491-7)</span> {#R1}

#### [3] <span>Butler, D.G., B.R. Cullis, A.R. Gilmour, and B.J. Gogel. 2009. **ASReml-R reference manual**. [Available Here](https://www.vsni.co.uk/downloads/asreml/release2/doc/asreml-R.pdf)</span> {#R2}

#### [4] <span>Endelman, J.B. 2011. Ridge Regression and Other Kernels for Genomic Selection with R Package rrBLUP. **The Plant Genome** 4: 250–255. [Available Here](http://dx.doi.org/10.3835/plantgenome2011.08.0024)</span> {#R3}

#### [5] <span>Pérez, P., and G. De Los Campos. 2014. Genome-wide regression & prediction with the BGLR statistical package. **Genetics** 198(2): 483–495.</span> {#R4}

#### [6] <span>Pérez-Elizalde, S., J. Cuevas, P. Pérez-Rodríguez, and J. Crossa. 2015. Selection of the Bandwidth Parameter in a Bayesian Kernel Regression Model for Genomic-Enabled Prediction. **Journal of Agricultural, Biological, and Environmental Statistics** 20(4): 512–532. [Available Here](http://link.springer.com/10.1007/s13253-015-0229-y)</span> {#R5}

#### [7] <span>Pritchard, J.K., M. Stephens, and P. Donnelly. 2000. Inference of Population Structure Using Multilocus Genotype Data. **Genetics** 155(2): 945-959. </span> {#R6}

#### [8] <span>Rutkoski, J.E., J. Poland, J.-L. Jannink, and M.E. Sorrells. 2013. Imputation of Unordered Markers and the Impact on Genomic Selection Accuracy. **G3: Genes, Genomes, Genetics** 3(3): 427–439. [Available Here](http://www.g3journal.org/content/3/3/427)</span> {#R7}

#### [9] <span>VanRaden, P.M. 2008. Efficient methods to compute genomic predictions. **Journal of dairy science** 91(11): 4414–23. [Available Here](http://www.sciencedirect.com/science/article/pii/S0022030208709901)</span> {#R8}

#### [10] <span>Yang, J., B. Benyamin, B.P. Mcevoy, S. Gordon, A.K. Henders, R. Dale, P.A. Madden, A.C. Heath, N.G. Martin, G.W. Montgomery, M.E. Goddard, and P.M. Visscher. 2010. Common SNPs explain a large proportion of heritability for human height. **Nature** 569: 565–569. [Available Here](http://www.nature.com/ng/journal/v42/n7/full/ng.608.html?foxtrotcallback=true)</span> {#R9}