  %\VignetteEngine{knitr::knitr}
  %\VignetteIndexEntry{How to use rehh}
  %\usepackage[utf8]{inputenc}
---
title: "Vignette for the package *rehh* (version 2+)"
author: "Mathieu Gautier, Alexander Klassmann and Renaud Vitalis"
date: "24/10/2016"
output:
  html_document:
    toc: yes
  pdf_document:
    fig_caption: yes
    number_sections: yes
    toc: yes
csl: genetics.csl
bibliography: vignette.bib
---

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

\newpage
This vignette is aimed at presenting additional information on the R package \emph{rehh} by describing how to use it to perform whole genome scan for footprints of selection using statistics related to the Extended Haplotype Homozygosity (EHH) [@Sabeti2002]. Importantly, the current implementation of tests assumes markers are bi-allelic.

The \emph{rehh} package is currently available for most platforms (Linux, MS Windows and MacOSX) from the CRAN repository (\url{http://cran.r-project.org/}) and may be installed using standard procedure. Once the package has been successfully installed on your system, it can be loaded using the following command:
```{r message=FALSE}
library(rehh)
```

#Input Files
\label{make.example.files}

The package \emph{rehh} basically requires as input:

1. haplotype data file(s) (in three possible format) for each population of interest (see \ref{hap.file})
2. a SNP information file (see \ref{map.file})

**Important Note: For a given chromosome, SNPs are assumed to be ordered in the same way in both the haplotype file (columns) and the SNP information file.**

For illustration purposes, example files that originate from a previously published study on the Creole cattle breed from Guadeloupe (CGU) [@Gautier2011] are provided in the package and can be copied in the working directory with the following command:
```{r results='hide'}
make.example.files()
```

In the following examples, the command \texttt{make.example.file()} is assumed to have been run (see above) so that example files are in the working directory.

##Haplotype data file
\label{hap.file}

Three haplotype input file formats are supported:

1.  a standard haplotype format. Each line represents a haplotype (the first element being an haplotype identifier) with SNP genotype in columns as in the example file \texttt{bta12\_cgu.hap} created when running the function \texttt{make.example.files()} and that consists of 280 haplotypes (identifier 1 to 280) of 1424 SNPs each. See \ref{LoadDataEx1} for a detailed example.
2.  A (transposed) format with haplotype in columns and SNPs in row as in the example file \texttt{bta12\_cgu.thap}  created when running the function \texttt{make.example.files()}. Note that this format is similar to the one produced by the phasing program \emph{SHAPEIT2} [@OConnell2014].  See \ref{LoadDataEx2} for a detailed example.
3.  the output file format from the phasing program \emph{fastPHASE} [@Scheet2006] as in the \texttt{bta12\_hapguess\_switch.out} example file  created when running the function \texttt{make.example.files()}. Note that haplotypes might originate from several different populations (i.e., if the -u \emph{fastPHASE} option was used).  See \ref{LoadDataEx3} for a detailed example.

By default alleles are assumed to be coded as 0 (missing data), 1 (ancestral allele) or 2 (derived allele). Recoding of the alleles in this format, according to the SNP information data file (see \ref{map.file}) can be performed with the \texttt{recode.allele} option of the function \texttt{data2haplohh()} (see \ref{LoadData}).

##SNP information data file
\label{map.file}

This data file should contain SNP information as in the \texttt{map.inp} example file created when running the function \texttt{make.example.files()}. Each line contains five columns corresponding to:

1. the SNP name
2. the SNP chromosome (or scaffold) of origin
3. the SNP position on the chromosome (or scaffold). Note that it is up to the user to choose either physical or genetic map positions (if available).
4. the SNP ancestral allele (as coded in the haplotype input file)
5. the SNP derived alleles (as coded in the haplotype input file)

The fourth and fifth columns (allele coding) should be filled in but the corresponding information is only used when activating the \texttt{recode.allele} option of the function \texttt{data2haplohh()} (see \ref{LoadData}). In that case, for each SNP, the allele specified in the fourth (respectively fifth) column will be recoded as 1 (respectively 2), any other allele name will be recoded as 0 (i.e., missing data). More importantly, it should be noticed that the ancestral or derived allele information associated to this coding are only relevant for within population tests (based on \emph{iHS}). In other words, if one is only interested in across-population tests (based on \emph{Rsb} or \emph{xpEHH}), assignment of the two SNP alleles in the fourth and fifth column may be performed randomly. 

As an illustration, the following R command displays the first five row of the \texttt{map.inp} example file created when running the function \texttt{make.example.files()}:

```{r}
head(read.table("map.inp"))
```


##Loading data files
\label{LoadData}
The \texttt{data2haplohh()} function allows to convert data file into an R object of class \texttt{haplohh} subsequently used by the functions of the \emph{rehh} package. The following main options are available to recode alleles or select SNPs (based on Minor Allele Frequency or percentage of missing data) and haplotypes (based on percentage of missing data): 

1. **Allele recoding option**: This option is activated with \texttt{recode.allele=TRUE} and allows to recode haplotype according to the ancestral and derived allele definition available in the SNP information file (fourth and fifth columns) as: 0 (missing data), 1 (ancestral allele) or 2 (derived allele). 
2. **Discard haplotypes with a high amount of missing data**. If haplotypes contain missing data (which is generally not the case since most phasing programs allow imputing missing genotypes), it is possible to discard those with less \texttt{min\_perc\_geno.hap} \% of the SNPs genotyped. By default \texttt{min\_perc\_geno.hap=100} meaning that only completely phased haplotypes are retained.   
3. **Discard SNPs with a high amount of missing data**. It is possible to discard SNPs genotyped on less than \texttt{min\_perc\_geno.snp} \% of the haplotypes. By default \texttt{min\_perc\_geno.snp=100} meaning that only fully genotyped SNPs are retained.
4. **Discard SNPs with a low Minor Allele Frequency**. It is possible to discard SNPs with a MAF below \texttt{min\_maf}. This is generally not recommended and by default \texttt{min\_maf=0} meaning that all SNPs are retained.\footnote{The arguments \texttt{min\_perc\_geno.hap}, \texttt{min\_perc\_geno.snp} and \texttt{min\_maf} are evaluated in this order.}

More details about the different arguments of the function are available in the documentation accessible using the command:
```{r}
?data2haplohh
```

In the following we detail three different examples based on the example data files provided with the package (see \ref{make.example.files}).

###Example 1: reading haplotype file in standard format
\label{LoadDataEx1}
In this example, the example haplotype input file \texttt{bta12\_cgu.hap} (standard format) and SNP information input files \texttt{map.inp} are converted into an \texttt{haplohh} object named \texttt{hap}. Because the map file contains information for SNPs mapping to other chromosomes than the one of interest (BTA12), we use the option \texttt{chr.name=12}. Allele recoding is activated (\texttt{recode.allele=TRUE}) to allow recoding alleles in the \emph{rehh} format (0,1 or 2).
```{r }
hap<-data2haplohh(hap_file="bta12_cgu.hap",map_file="map.inp",
                  recode.allele=TRUE,chr.name=12)
```

If no value is specified for \texttt{chr.name} argument and more than one chromosome is detected in the map file, the function asks interactively which chromosome to choose:

```{r eval=FALSE}
hap<-data2haplohh(hap_file="bta12_cgu.hap",map_file="map.inp",
                  recode.allele=TRUE) 
```
```{r echo=FALSE}
cat("More than one chromosome name in Map file: map.inp\nWhich chromosome should be considered among:\n1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29\n1:")
```
```{r eval=FALSE}
12
```
```{r echo=FALSE}
hap<-data2haplohh(hap_file="bta12_cgu.hap",map_file="map.inp",recode.allele=TRUE,chr.name=12)
```

Finally, as an example of error message, the following message is prompted if the number of SNPs in the chromosome (for instance when a wrong chromosome is declared) does not correspond to the one in the haplotype file:

```{r error=TRUE}
hap<-data2haplohh(hap_file="bta12_cgu.hap",map_file="map.inp",
                  recode.allele=TRUE,chr.name=18)
```

###Example 2: reading haplotype file in transposed format (\emph{SHAPIT2}--like)
\label{LoadDataEx2}

In this example, the example haplotype input file \texttt{bta12\_cgu.thap} (transposed format) and SNP information input files \texttt{map.inp} are converted into an \texttt{haplohh} object named \texttt{hap}. Setting \texttt{haplotype.in.columns=TRUE} informs the function that the haplotype file is in transposed format:
```{r}
hap<-data2haplohh(hap_file="bta12_cgu.thap",map_file="map.inp",haplotype.in.columns=TRUE,
                  recode.allele=TRUE,chr.name=12)
```

###Example 3: reading haplotype file in \emph{fastPHASE} output format
\label{LoadDataEx3}
In this example, the example \emph{fastPHASE} output file \texttt{bta12\_hapguess\_switch.out} and SNP information input files \texttt{map.inp} are 
converted into a \texttt{haplohh} object named \texttt{haplo}. As explained above we use the option \texttt{chr.name=12}. 
Because, haplotypes originate from several populations (the -u \emph{fastPHASE} option was used), we specify the population of interest 
(in our example the 280 haplotypes from the CGU population, see above) using the option \texttt{popsel=7} 
(7 corresponding to the code of CGU in the example \emph{fastPHASE} input files).

```{r}
hap<-data2haplohh(hap_file="bta12_hapguess_switch.out",map_file="map.inp",
                  recode.allele=TRUE,popsel=7,chr.name=12)
```

If no value is specified for the \texttt{popsel} argument and more than one population is detected in the \emph{fastPHASE} output file, the function asks interactively which population to chose:

```{r eval=FALSE}
hap<-data2haplohh(hap_file="bta12_hapguess_switch.out",map_file="map.inp",
                  recode.allele=TRUE,chr.name=12)
```
```{r echo=FALSE}
cat("Map file seems OK: 1424  SNPs declared for chromosome 12\nLooks like a FastPHASE haplotype file\nHaplotypes originate from  8  different populations in the fastPhase output file\nChosen pop. is not in the list of pop. number: 1 2 3 4 5 6 7 8\nWhich population should be considered among: 1 2 3 4 5 6 7 8\n1: 
")
```
```{r eval=FALSE}
7
```
```{r echo=FALSE}
hap<-data2haplohh(hap_file="bta12_hapguess_switch.out",map_file="map.inp",recode.allele=TRUE,popsel=7,chr.name=12)
```

#Computing \emph{EHH}-based statistics for individual markers with the \texttt{calc\_ehh()}, \texttt{calc\_ehhs()} and \texttt{scan\_hh()} functions
##Definition and Computation

###The (allele-specific) Extended Haplotype Homozygosity (\emph{EHH})
\label{EHH}
For a given core allele (i.e., the ancestral or derived allele) at a focal SNP, the (allele--specific) extended haplotype homozygosity (\emph{EHH}) is defined as the probability that two randomly chosen chromosomes (carrying the core allele considered) are identical by descent (as assayed by homozygosity at all SNPs) over a given surrounding chromosome region [@Sabeti2002]. The \emph{EHH} thus aims at measuring to which extent an extended haplotype is transmitted without recombination. In practice, the \emph{EHH} ($\mathrm{EHH}_{s,t}$) of a tested core allele $a_s$ ($a_s=1$ or $a_s=2$) for a focal SNP $s$ over the chromosome interval extending to SNP $t$ is computed as:
\begin{equation}
\label{ehh.def}
\mathrm{EHH}_{s,t}=\frac{1}{n_{a_s}(n_{a_s}-1)}\sum\limits_{k=1}^{K_{a_s,t}}n_k(n_k-1)
\end{equation}
where $n_{a_s}$ represents the number of haplotype carrying the core allele $a_s$, $K_{a_s,t}$ represents the number of different extended haplotypes (from SNP $s$ to SNP $t$) carrying $a_s$ and $n_k$ is the number of the extended haplotype $k$ (i.e., $n_{a_s}=\sum\limits_{k=1}^{K_{a_s,t}}n_k$).

###The integrated (allele-specific) \emph{EHH} (\emph{iHH})
\label{iHH}
By definition, irrespective of the allele considered, \emph{EHH} starts at 1, and decays monotonically to 0 with increasing distance from the focal SNP. For a given core allele, the integrated \emph{EHH} (\emph{iHH}) is defined as the area under the \emph{EHH} curve with respect to map position [@Voight2006]\footnote{In their seminal paper, Voight et al. considered genetic distances and apply a penalty (proportional to physical distances) for successive SNPs separated by more than 20 kb. In addition, they did not compute \emph{iHH} if any physical distance between a pair of neighboring SNPs was above 200 kb. We did not implement such an approach in \emph{rehh} although this might easily be done (providing relevant information is available for the genome of interest) by modifying the positions of the markers in SNP information input file.}. In \emph{rehh}, \emph{iHH} is computed using the trapezoid method. In practice, the integral might only be computed for the regions of the curve above an arbitrarily small \emph{EHH} value (e.g., \emph{EHH}>0.05). 

###The site-specific Extended Haplotype Homozygosity (\emph{EHHS})
\label{EHHS}
For a given core SNP, the (site--specific) extended haplotype homozygosity (\emph{EHHS}) is defined as the probability that two randomly chosen chromosomes are identical by descent (as assayed by homozygosity at all SNPs) over a given surrounding chromosome region. \emph{EHHS}  might roughly be viewed as linear combination of the \emph{EHH}'s for the two alternative alleles with weights function of the corresponding allele frequencies. Two different \emph{EHHS} estimators further denoted as $\mathrm{EHHS}^{\text{Sab}}$ [@Sabeti2007] and $\mathrm{EHHS}^{\text{Tang}}$ [@Tang2007] have been proposed. For a focal SNP $s$ over a chromosome interval extending to SNP $t$, these are computed as (following the same notation as above):
\begin{equation}
\label{ehhssab.def}
\mathrm{EHHS}^{\text{Sab}}_{s,t}=\frac{1}{n_s(n_s-1)}\sum\limits_{a_s=1}^{a_s=2}\left(\sum\limits_{k=1}^{K_{a_s,t}}n_k(n_k-1)\right)
\end{equation}
where $n_s=\sum\limits_{a_s=1}^{a_s=2}n_{a_s}$ and
\begin{equation}
\label{ehhstang.def}
\mathrm{EHHS}^{\text{Tang}}_{s,t}=\frac{1-h_{hap}^{(s,t)}}{1-h_{all}^{(s)}}
\end{equation}
where:

1. $h_{all}^{(s)}=\frac{n_s}{n_s-1}\left(1-\frac{1}{n_s^2}\sum\limits_{a_s=1}^{a_s=2}n_{a_s}^2 \right )$ is an estimator of the focal SNP heterozygosity
2. $h_{hap}^{(s,t)}=\frac{n_s}{n_s-1}\left(1-\frac{1}{n_s^2}\sum\limits_{a_s=1}^{a_s=2}\left(\sum\limits_{k=1}^{K_{a_s,t}}n_k^2 \right )\right)$ is an estimator of haplotype heterozygosity across the chromosome region extending from SNP $s$ to SNP $t$. 

###The integrated \emph{EHHS} (\emph{iES})
\label{iES}
As for the \emph{EHH} (see \ref{iHH}), \emph{EHHS} starts at 1 and decays monotonically to 0 with increasing distance from the focal SNP. For a given focal SNP, and in a similar fashion as the \emph{iHH}, \emph{iES} is defined as the integrated \emph{EHHS} [@Tang2007]. Depending on the \emph{EHHS} estimator considered (respectively, $\mathrm{EHHS}^{\text{Sab}}$ and $\mathrm{EHHS}^{\text{Tang}}$), two different \emph{iES} estimators, that we further denoted as $\mathrm{iES}^{\text{Sab}}$ and $\mathrm{iES}^{\text{Tang}}$ respectively, can be computed. As for \emph{iHH}, the  \emph{iES} integral is computed using the trapezoid method and might only be computed for the regions of the curve above an arbitrarily small \emph{EHHS} value (e.g., \emph{EHHS}>0.05). 

###Dealing with missing data

In the computation of both \emph{EHH} and \emph{EHHS} from a focal SNP $s$ to a SNP $t$, only extended haplotypes with no missing data are considered. As a consequence, the number of extended haplotypes retained to compute these two statistics might decrease with distance from the focal SNP. However if the number of available extended haplotypes falls below a threshold (e.g., \texttt{limhaplo}=5), \emph{EHH} and \emph{EHHS} are not computed further. Note however that most phasing programs (such as \emph{fastPHASE} or \emph{SHAPEIT2}) allow to impute missing genotypes resulting in phased haplotypes with no missing data.

##The function \texttt{calc\_ehh()}
\label{calcehh}
The \texttt{calc\_ehh()} function allows to compute \emph{EHH} for both the ancestral ($a_s=1$) and derived ($a_s=2$) alleles at the $s^\text{th}$ SNP relative to each SNP ($t$) upstream and downstream and corresponding \emph{iHH}. The two options \texttt{limehh} and \texttt{limhaplo} allow to specify condition to stop computing \emph{EHH} (see \ref{EHH}). By default \texttt{limehh}=0.05 and \texttt{limhaplo}=2. Finally, if \texttt{plotehh=TRUE}, the decay of \emph{EHH} for both the ancestral and derived alleles is plotted against SNP map position (\texttt{main\_leg} allows to change the plot legend). More details are available in the R documentation by using the command:
```{r eval=FALSE}
?calc_ehh
```
In the following example,the \emph{EHH} statistics are computed for both the ancestral and derived allele of the $456^{\text{th}}$ focal SNP. Note that the \texttt{haplohh\_cgu\_bta12} object was generated using the \texttt{data2haplohh()} function with the example input files (\ref{LoadDataEx1}). For convenience, it is stored as an example object (accessible with the R function \texttt{data}) as shown below:
```{r, fig.align='center',fig.keep='all'}
#example haplohh object (280 haplotypes, 1424 SNPs) see ?haplohh_cgu_bta12 for details
data(haplohh_cgu_bta12)
#computing EHH statistics for the focal SNP at position 456 
#which display a strong signal of selection
res.ehh<-calc_ehh(haplohh_cgu_bta12,mrk=456) 
```
The five different elements of the resulting \texttt{res.ehh} object are as follows:
```{r}
res.ehh$ehh[1:2,454:458]
```
```{r}
res.ehh$nhaplo_eval[1:2,454:458]
```
```{r}
res.ehh$freq_all1
```
```{r}
res.ehh$ihh
```
In addition, as \texttt{plotehh=TRUE} by default, we obtain the following plot (Figure \ref{fig:ehhplot}):
```{r ehhplot,echo=FALSE,fig.align='center',,out.height="7cm",fig.cap='Graphical output for the \\texttt{calc\\_ehh()} function',fig.pos='!h',fig.lp='fig:'}
sel_reg <- (colSums(res.ehh$nhaplo_eval) > 0) ; mrk=456
matplot(haplohh_cgu_bta12@position[sel_reg]/1000000, t(res.ehh$ehh[, sel_reg]),col = c("red", "blue"), lty = 1, type = "l",main = "", bty = "n", xlab = "Position (Mb)",ylab = expression(Extended ~ haplotype ~ homozygosity ~ (italic(EHH))))
abline(v = haplohh_cgu_bta12@position[mrk]/1000000, lty = 2)
legend("topleft", c("Ancestral Allele", "Derived Allele"),col = c("red", "blue"), bty = "n", lty = 1)
```

##The function \texttt{calc\_ehhs()}

The \texttt{calc\_ehhs()} function allows to compute the \emph{EHHS} (both the $\mathrm{EHHS}^{\text{Sab}}$ and $\mathrm{EHHS}^{\text{Tang}}$ estimators) at the $s^\text{th}$ SNP relative to each SNP ($t$) upstream and downstream. This function also compute the corresponding \emph{iES} ($\mathrm{iES}^{\text{Sab}}$ and $\mathrm{iES}^{\text{Tang}}$ estimators respectively). The two options \texttt{limehhs} and \texttt{limhaplo} allow to specify condition to stop computing \emph{EHHS} (see \ref{EHHS}). By default \texttt{limehhs}=0.05 and \texttt{limhaplo}=2. Finally, if \texttt{plotehhs=TRUE}, the decay of \emph{EHHS} is plotted against SNP map position (\texttt{main\_leg} allows to change the plot legend). More details are available in the R documentation by using the command:
```{r, eval=FALSE}
?calc_ehhs
```
In the following example, the \emph{EHHS} statistics are computed for the $456^{\text{th}}$ focal SNP on the \texttt{haplohh\_cgu\_bta12} object defined above (see \ref{calcehh}) was generated using the \texttt{data2haplohh()} function with the example input files (see \ref{LoadDataEx1}) described above. For convenience, it is stored as an example object (accessible with the R function \texttt{data}).

```{r}
#example haplohh object (280 haplotypes, 1424 SNPs) see ?haplohh_cgu_bta12 for details
data(haplohh_cgu_bta12)
#computing EHH statistics for the focal SNP at position 456
#which display a strong signal of selection
res.ehhs<-calc_ehhs(haplohh_cgu_bta12,mrk=456)
```
The five different elements of the resulting \texttt{res.ehhs} object are as follows:
```{r}
res.ehhs$EHHS_Sabeti_et_al_2007[453:459] 
```
```{r}
res.ehhs$EHHS_Tang_et_al_2007[453:459] 
```
```{r}
res.ehhs$nhaplo_eval[453:459] 
```
```{r}
res.ehhs$IES_Tang_et_al_2007
```
```{r}
res.ehhs$IES_Sabeti_et_al_2007
```
In addition, as \texttt{plotehh=TRUE} by default, we obtain the following plot (Figure \ref{fig:ehhsplot}):
```{r ehhsplot,echo=FALSE,fig.align='center',out.height="7cm",fig.cap='Graphical output for the \\texttt{calc\\_ehhs()} function',fig.pos='!h',fig.lp='fig:'}
sel_reg <- (colSums(res.ehh$nhaplo_eval) > 0) ; mrk=456
plot(haplohh_cgu_bta12@position[sel_reg]/1000000, res.ehhs$EHHS_Sabeti_et_al_2007[sel_reg],col = "red", lty = 1, type = "l",main = "", bty = "n", xlab = "Position (Mb)",ylab = expression(Site ~ specific ~ italic(EHH) ~ (italic(EHHS)) ))
lines(haplohh_cgu_bta12@position[sel_reg]/1000000, res.ehhs$EHHS_Tang_et_al_2007[sel_reg], col = "blue")
abline(v = haplohh_cgu_bta12@position[mrk]/1000000, lty = 2)
legend("topleft", c("Sabeti et al. (2007)", "Tang et al. (2007)"),col = c("red", "blue"), bty = "n", lty = 1)
```

##The function \texttt{scan\_hh()}
\label{fun.scanhh}
The \texttt{scan\_hh()} function allows to efficiently compute \emph{IHH} (for both the ancestral and derived alleles) and \emph{IES} (both the $\mathrm{iES}^{\text{Sab}}$ and $\mathrm{iES}^{\text{Tang}}$ estimators) for all the SNPs in the \texttt{haplohh} object considered. The options \texttt{limehh}, \texttt{limehhs} and \texttt{limhaplo} specify conditions to stop computing \emph{EHH} and \emph{EHHS}. By default \texttt{limehh}=\texttt{limehhs}=0.05 and \texttt{limhaplo}=2. Finally, the option \texttt{threads}, set by dafault to \texttt{threads}=1, allows to specify the number of available threads to parallelize computation (parallelization being carried out over SNPs).
For instance to scan the \texttt{haplohh\_cgu\_bta12} object (containing data on 1424 SNPs for 280 haplotypes), one may use the following command:
```{r,eval=TRUE}
data(haplohh_cgu_bta12)
res.scan<-scan_hh(haplohh_cgu_bta12)
```
The resulting object \texttt{res.scan} is a data frame with \texttt{haplohh\_cgu\_bta12\@nsnp} (number of SNPs declared in the \texttt{haplohh} object) and seven columns giving for each SNP:

1. The SNP chromosome of origin
2. The SNP position
3. The SNP ancestral allele frequency
4. The estimated \emph{IHH} for the ancestral allele (\emph{IHH\_A})
5. The estimated \emph{IHH} for the derived allele (\emph{IHH\_D})
6. The estimated $\mathrm{iES}^{\text{Tang}}$
7. The estimated $\mathrm{iES}^{\text{Sab}}$.

As an example, the following R codes provide the dimension and the first five rows of the \texttt{res.scan} data frame obtained above:
```{r}
dim(res.scan)
```
```{r}
head(res.scan)
```
Note that running \texttt{scan\_hh()} is more efficient than running \texttt{calc\_ehh()} and \texttt{calc\_ehhs()} in turn as shown in the example below (\texttt{scan\_hh()}:.
```{r}
system.time(res.scan<-scan_hh(haplohh_cgu_bta12))
```
```{r}
foo<-function(haplo){
res.ihh=res.ies=matrix(0,haplo@nsnp,2)
for(i in 1:length(haplo@position)){
  res.ihh[i,]=calc_ehh(haplo,mrk=i,plotehh=FALSE)$ihh
  tmp=calc_ehhs(haplo,mrk=i,plotehhs=FALSE)
  res.ies[i,1]=tmp$IES_Tang_et_al_2007
  res.ies[i,2]=tmp$IES_Sabeti_et_al_2007  
}
list(res.ies=res.ies,res.ihh=res.ihh)
}
system.time(res.scan2<-foo(haplohh_cgu_bta12))
```
Note however that the same results are obtained (since the same options were used) as illustrated by the following R code:
```{r}
sum(res.scan2$res.ihh[,1]!=res.scan[,4]) + sum(res.scan2$res.ihh[,2]!=res.scan[,5]) +
sum(res.scan2$res.ies[,1]!=res.scan[,6]) + sum(res.scan2$res.ies[,2]!=res.scan[,7])
```
#Computing \emph{iHS}, \emph{Rsb} and \emph{xpEHH}: the \texttt{ihh2ihs()},\texttt{ies2rsb()} and \texttt{ies2xpehh()} functions

##Within population test: the \emph{iHS}
###Definition
\label{def.ihs}
Let $\mathrm{UniHS}$ represent the log-ratio of the  \emph{iHH} for its ancestral (\emph{iHH}$_a$) and derived (\emph{iHH}$_d$) alleles:
$$\mathrm{UniHS}=\log\left(\frac{\mathrm{iHH}_a}{\mathrm{iHH}_d}\right)$$
The \emph{iHS} of a given focal SNP $s$ ($\mathrm{iHS}(s)$) is then defined as its standardized $\mathrm{UniHS}$ ($\mathrm{UniHS}(s)$) following [@Voight2006]:
$$\mathrm{iHS}(s)=\frac{\mathrm{UniHS}(s) - \mu^{p_s}_\mathrm{UniHS}}{\sigma^{p_s}_\mathrm{UniHS}}$$
where $\mu^{p_s}_\mathrm{UniHS}$ and $\sigma^{p_s}_\mathrm{UniHS}$ represent the average and standard deviation of the  $\mathrm{UniHS}$ computed over all the SNPs with a derived allele frequency $p_s$ similar to that of the core SNP $s$. In practice, the derived allele frequencies are generally binned so that each bin are large enough (e.g., >10 SNPs) to obtain reliable estimate of  $\mu^{p_s}_\mathrm{UniHS}$ and $\sigma^{p_s}_\mathrm{UniHS}$.

Note that the \emph{iHS} is constructed to have an approximately standard Gaussian distribution and to be comparable across SNPs regardless of their underlying allele frequencies. Hence, one may further transform \emph{iHS} into $p_\mathrm{iHS}$ [@Gautier2011]:
$$p_\mathrm{iHS}=-\log_{10}\left(1-2|\Phi\left(\mathrm{iHS}\right)-0.5|\right)$$
where $\Phi\left(x\right)$ represents the Gaussian cumulative distribution function. Assuming most of the genotyped SNPs behave neutrally (i.e., the genome-wide empirical \emph{iHS} distribution is a fair approximation of the neutral distribution), $p_\mathrm{iHS}$ might thus be interpreted as a two-sided P-value (on a $-\log_{10}$ scale) associated to the neutral hypothesis of no selection.

###The function \texttt{ihh2ihs()}
\label{cgu.ihs}
The \texttt{ihh2ihs()} function allows to compute \emph{iHS} using a matrix of \emph{iHH} statistics (for both the ancestral and derived alleles) in the same format as obtained from the \texttt{scan\_hh()} function (see \ref{fun.scanhh}). The argument \texttt{minmaf} allows to remove SNPs according to their MAF (by default SNPs with a MAF<\texttt{minmaf}=0.05 are discarded from the standardization). The argument \texttt{freqbin} controls the size of the allele frequency bins used to perform standardization (see \ref{def.ihs}). More precisely allele frequency bins vary from \texttt{minmaf} to 1-\texttt{minmaf} per step of size \texttt{freqbin} (by default \texttt{freqbin}=0.025).  Note that if \texttt{freqbin} is set to 0 (e.g., with a large number of SNPs and few haplotypes), standardization is performed considering each observed frequency as a frequency class.

For instance, to perform a whole genome scan one might run \texttt{scan\_hh()} in turn on haplotype data from each chromosome and concatenate the resulting matrices before standardization. In the following example, we assume that the haplotype files are named as \texttt{hap\_chr\_i.pop1} where the chromosome number $i$ goes from 1 to 29 and the SNP information file is named \texttt{snp.info}. The R code below then generates a matrix \texttt{wg.res} with \emph{iHH$_a$} and \emph{iHH$_d$} estimates for all SNPs in an appropriate format to perform standardization with the \texttt{ihh2ihs} function:
```{r,eval=FALSE}
for(i in 1:29){
 hap_file=paste("hap_chr_",i,".pop1",sep="")
 data<-data2haplohh(hap_file="hap_file","snp.info",chr.name=i)
 res<-scan_hh(data)
 if(i==1){wg.res<-res}else{wg.res<-rbind(wg.res,res)}
   }
wg.ihs<-ihh2ihs(wg.res)
```
As a matter of illustration, results of a similar genome scan [@Gautier2011] are provided as example data sets. The following R code allows to compute the \emph{iHS} for the CGU population:
```{r}
data(wgscan.cgu)
## results from a genome scan (44,057 SNPs) see ?wgscan.eut and ?wgscan.cgu for details
ihs.cgu<-ihh2ihs(wgscan.cgu)
```
The corresponding object \texttt{ihs.cgu} is a list with two elements corresponding to

1. a matrix of SNP \emph{iHS} and the corresponding $p_\mathrm{iHS}$  (P-values in a $-\log_10$ scale assuming the \emph{iHS} are normally distributed under the neutral hypothesis).  For instance, the five first rows of the \texttt{ihs.cgu\$iHS} data frame are displayed below using the following R command:
```{r}
head(ihs.cgu$iHS)
```
2.a matrix summarizing the allele frequency bins. For instance, the five first rows of the \texttt{ihs.cgu\$frequency.class} data frame are displayed below using the following R command:
```{r}
head(ihs.cgu$frequency.class)
```

###Manhattan plot of the results: the function \texttt{ihsplot()}

The \texttt{ihsplot()} function allows to draw a Manhattan plot of the Whole Genome scan results as stored in the list object produced by the function \texttt{ihh2ies()}.  Various options are available to modify the graphical display (see \texttt{?ihsplot}).
```{r,eval=FALSE}
ihsplot(ihs.cgu,plot.pval=TRUE,ylim.scan=2,main="iHS (CGU cattle breed)") 

```

```{r ,echo=FALSE,fig.ext="png",fig.align='center',fig.width=16,fig.height=12,fig.cap='Graphical output for the \\texttt{ihsplot()} function'}
layout(matrix(1:2,2,1))
plot.pval = TRUE; ylim.scan = 2; pch = 16; cex = 0.75; cex.lab = 1; main = "iHS (CGU cattle breed)"; cex.main = 1.;cex.axis=0.7
data = ihs.cgu$iHS
lst_chrm <- unique(data$CHR) ; nbr_chrm <- length(lst_chrm) ; cum <- rep(0, (nbr_chrm + 1))
if (nbr_chrm > 1) {
  for (i in 1:nbr_chrm) {cum[i + 1] <- tail(data$POSITION[data$CHR == lst_chrm[i]],n = 1)}
}
cum <- cumsum(cum) ; pos <- rep(0, length(data$CHR)) ; pos_labels <- rep(0, nbr_chrm)
for (i in 1:nbr_chrm) {
  pos[data$CHR == lst_chrm[i]] <- data$POSITION[data$CHR ==lst_chrm[i]] + cum[which(lst_chrm == lst_chrm[i])]
  pos_labels[i] <- (cum[i] + cum[i + 1])/2
}
col_chr <- data$CHR
par(mar = c(5, 5, 4, 2) + 0.1)
plot(pos, data$iHS, pch = pch, cex = cex, las = 1, col = col_chr,xaxt = "n", xlab = "Chromosome", ylab = expression(italic(iHS)), 
     cex.lab = cex.lab, main = main, cex.main = cex.main,cex.axis=cex.axis)
axis(1, at = pos_labels, labels = lst_chrm, las = 1,cex.lab=cex.lab,cex.axis=cex.axis)
abline(h = c(-ylim.scan, ylim.scan), lty = 2)
par(mar = c(5, 5, 4, 2) + 0.1)
plot(pos, data$"-log10(p-value)", pch = pch, cex = cex,las = 1, col = col_chr, xaxt = "n", xlab = "Chromosome", 
         ylab = expression("-" * log[10] * "[1" ~ "-" ~"2" ~ "|" ~ Phi[scriptstyle(italic(iHS))] ~"-" ~ 0.5 * "|]"),            cex.lab = cex.lab, main = main, cex.main = cex.main,cex.axis=cex.axis)
axis(1, at = pos_labels, labels = lst_chrm, las = 1,cex.lab=cex.lab,cex.axis=cex.axis)
abline(h = c(-ylim.scan, ylim.scan), lty = 2)
```

##The \emph{Rsb}--based pairwise population test 
###Definition
For a given SNP $s$, let $$\mathrm{LRiES}(s)^{\text{Tang}}=\log\left(\frac{\mathrm{iES}_\text{pop1}(s)^{\text{Tang}}}{\mathrm{iES}_\text{pop2}(s)^{\text{Tang}}}\right)$$ represent the log-ratio of the $\mathrm{iES}_\text{pop1}(s)^{\text{Tang}}$ and $\mathrm{iES}_\text{pop2}(s)^{\text{Tang}}$ computed in the pop1 and pop2 populations (see \ref{iES}).

The \emph{Rsb} for a given focal SNP is then defined as the standardized  $\mathrm{LRiES}(s)^{\text{Tang}}$ [@Tang2007]:

\begin{equation}
\mathrm{rSB}(s)=\frac{\mathrm{LRiES}(s)^{\text{Tang}} - \text{med}_{\mathrm{LRiES}^{\text{Tang}}}}{\sigma_{\mathrm{LRiES}^{\text{Tang}}}}
\end{equation}
where $\text{med}_{\mathrm{LRiES}^{\text{Tang}}}$ and $\sigma_{\mathrm{LRiES}^{\text{Tang}}}$ represent the median and standard deviation of the $\mathrm{LRiES}(s)^{\text{Tang}}$ computed over all the analyzed SNPs. Note that the median is used instead of the mean because it is less sensitive to extreme data points [@Tang2007]. 
More  importantly, it should be noticed that the information about the ancestral and derived status of alleles at the focal SNP is not needed.

As for the \emph{iHS} (see \ref{def.ihs}), \emph{Rsb} is constructed to have an approximately standard Gaussian distribution and may further be transformed into $p_\mathrm{rSB}$:
\begin{equation}
p_\mathrm{rSB}=-\log_{10}\left(1-2|\Phi\left(\mathrm{rSB}\right)-0.5|\right)
\end{equation}
where $\Phi\left(x\right)$ represents the Gaussian cumulative distribution function. Assuming most of the genotyped SNPs behave neutrally (i.e., the genome-wide empirical \emph{Rsb} distribution is a fair approximation of their corresponding neutral distributions), $p_\mathrm{rSB}$ might thus be interpreted as a two-sided P-value (in a $-\log_{10}$ scale) associated to the neutral hypothesis of no selection.
Alternatively, $p_\mathrm{rSB}$ might also be computed [@Gautier2011]:
\begin{equation}
p\prime_\mathrm{rSB}=-\log_{10}\left(|\Phi\left(\mathrm{rSB}\right)|\right)
\end{equation}
$p\prime_\mathrm{rSB}$ and $p\prime_\mathrm{rSB}$ might then be interpreted as a one-sided P-value (in a $-\log_{10}$ scale) allowing the identification of those sites displaying a significantly high extended haplotype homozygosity in population $pop2$ (represented in the denominator of the corresponding $\mathrm{LRiES}$) relatively to the $pop1$ reference population.

###The function \texttt{ies2rsb()}
\label{res.rsb}

The \texttt{ies2rsb()} function allows to compute \emph{Rsb} using two data frames containing the \emph{iES} statistics for each of the two populations considered in the same format as the one obtained after running the \texttt{scan\_hh()} function (see \ref{fun.scanhh}).
For instance, to perform a genome scan one might first run for each population \texttt{scan\_hh()} in turn on haplotype data from each chromosome and concatenate the resulting matrices. In the following example, we assume that the haplotype files are named as \texttt{hap\_chr\_i.pop1} and  \texttt{hap\_chr\_i.pop2} where $i$ is the chromosome number (going from 1 to 29), the suffixes pop1 and pop2 indicate the population of origin and the SNP information file is named \texttt{snp.info}. The R code below then generates two data frames (\texttt{wg.res.pop1} and \texttt{wg.res.pop2}) containing the results from all SNPs in the appropriate format to compute \emph{Rsb} with the \texttt{ies2rsb()} function:

```{r,eval=FALSE}
for(i in 1:29){
 hap_file=paste("hap_chr_",i,".pop1",sep="")
 data<-data2haplohh(hap_file="hap_file","snp.info",chr.name=i)
 res<-scan_hh(data)
 if(i==1){wg.res.pop1<-res}else{wg.res.pop1<-rbind(wg.res.pop1,res)}
 hap_file=paste("hap_chr_",i,".pop2",sep="")
 data<-data2haplohh(hap_file="hap_file","snp.info",chr.name=i)
 res<-scan_hh(data)
 if(i==1){wg.res.pop2<-res}else{wg.res.pop2<-rbind(wg.res.pop2,res)}
 }
wg.rsb<-ies2rsb(wg.res.pop1,wg.res.pop2)
```
As a matter of illustration, one may consider results from a similar genome scan [@Gautier2011] provided as example data sets and compute for each SNP the \emph{Rsb} between the CGU and EUT populations as follows:
```{r}
data(wgscan.cgu) ; data(wgscan.eut)
## results from a genome scan (44,057 SNPs) see ?wgscan.eut and ?wgscan.cgu for details
cguVSeut.rsb<-ies2rsb(wgscan.cgu,wgscan.eut,"CGU","EUT")
```

The resulting object \texttt{cguVSeut.rsb} is a data frame with of SNP \emph{Rsb} (and corresponding P-Values assuming \emph{Rsb} are normally distributed under the neutral hypothesis). Note that either bilateral (default) or unilateral might be performed (\texttt{method} argument). The five first rows of the \texttt{cguVSeut.rsb} data frame are displayed below using the following R command:

```{r}
head(cguVSeut.rsb)
```
###Manhattan plot of the results: the function \texttt{rsbplot()}

The \texttt{rsbplot()} function allows to draw a Manhattan plot of the Whole Genome scan results as stored in the data frame produced by the function \texttt{ies2rsb()}.  Various options are available to modify the graphical display (see \texttt{?rsbplot}). As an example, the Figure \ref{fig:rsbplot} below provides the output of the function \texttt{rsbplot} for the \emph{xpEHH} computed above across the CGU and EUT populations (see \ref{res.rsb}).  Figure \ref{fig:rsbplot} was drawn using the following R code:
```{r,eval=FALSE}
rsbplot(cguVSeut.rsb,plot.pval=TRUE)
```

```{r rsbplot,echo=FALSE,fig.ext="png",fig.align='center',fig.width=16,fig.height=12,fig.lp='fig:',fig.cap='Graphical output for the \\texttt{rsbplot()} function',fig.pos='!h'}
layout(matrix(1:2,2,1))
plot.pval = TRUE; ylim.scan = 2; pch = 16; cex = 0.75; cex.lab = 1; main = "Rsb (CGU vs. EUT)"; cex.main = 1.;cex.axis=0.7
data = cguVSeut.rsb
lst_chrm <- unique(data$CHR) ; nbr_chrm <- length(lst_chrm) ; cum <- rep(0, (nbr_chrm + 1))
if (nbr_chrm > 1) {
  for (i in 1:nbr_chrm) {cum[i + 1] <- tail(data$POSITION[data$CHR == lst_chrm[i]],n = 1)}
}
cum <- cumsum(cum) ; pos <- rep(0, length(data$CHR)) ; pos_labels <- rep(0, nbr_chrm)
for (i in 1:nbr_chrm) {
  pos[data$CHR == lst_chrm[i]] <- data$POSITION[data$CHR ==lst_chrm[i]] + cum[which(lst_chrm == lst_chrm[i])]
  pos_labels[i] <- (cum[i] + cum[i + 1])/2
}
col_chr <- data$CHR
par(mar = c(5, 5, 4, 2) + 0.1)
plot(pos, data[,3], pch = pch, cex = cex, las = 1, col = col_chr,xaxt = "n", xlab = "Chromosome", ylab = expression(italic(Rsb)), 
     cex.lab = cex.lab, main = main, cex.main = cex.main,cex.axis=cex.axis)
axis(1, at = pos_labels, labels = lst_chrm, las = 1,cex.lab=cex.lab,cex.axis=cex.axis)
abline(h = c(-ylim.scan, ylim.scan), lty = 2)
par(mar = c(5, 5, 4, 2) + 0.1)
plot(pos, data$"-log10(p-value)", pch = pch, cex = cex,las = 1, col = col_chr, xaxt = "n", xlab = "Chromosome", 
         ylab = expression("-" * log[10] * "[1" ~ "-" ~"2" ~ "|" ~ Phi[scriptstyle(italic(Rsb))] ~"-" ~ 0.5 * "|]"), cex.lab = cex.lab, main = main, cex.main = cex.main,cex.axis=cex.axis)
axis(1, at = pos_labels, labels = lst_chrm, las = 1,cex.lab=cex.lab,cex.axis=cex.axis)
abline(h = c(-ylim.scan, ylim.scan), lty = 2)
```

##The \emph{xpEHH}--based pairwise population test 
###Definition
The \emph{xpEHH} statistics [@Sabeti2007] is similar to the \emph{Rsb} except that it is based on the $\mathrm{iES}_\text{pop2}(s)^{\text{Sab}}$ (instead of $\mathrm{iES}_\text{pop2}(s)^{\text{Tang}}$) estimator of the \emph{iES} (see \ref{iES}). Hence, for or a given SNP $s$, let $$\mathrm{LRiES}(s)^{\text{Sab}}=\log\left(\frac{\mathrm{iES}_\text{pop1}(s)^{\text{Sab}}}{\mathrm{iES}_\text{pop2}(s)^{\text{Sab}}}\right)$$ represent the log-ratio of the $\mathrm{iES}_\text{pop1}(s)^{\text{Sab}}$ and $\mathrm{iES}_\text{pop2}(s)^{\text{Sab}}$ computed in the pop1 and pop2 populations (see \ref{iES}).

The \emph{xpEHH} for a given focal SNP is then defined as the standardized  $\mathrm{LRiES}(s)^{\text{Sab}}$ [@Sabeti2007]:

\begin{equation}
\mathrm{rSB}(s)=\frac{\mathrm{LRiES}(s)^{\text{Sab}} - \text{med}_{\mathrm{LRiES}^{\text{Sab}}}}{\sigma_{\mathrm{LRiES}^{\text{Sab}}}}
\end{equation}
where $\text{med}_{\mathrm{LRiES}^{\text{Sab}}}$ and $\sigma_{\mathrm{LRiES}^{\text{Sab}}}$ represent the median and standard deviation of the $\mathrm{LRiES}(s)^{\text{Sab}}$ computed over all the analyzed SNPs. 
More  importantly, it should be noticed that the information about the ancestral and derived status of alleles at the focal SNP is not needed.

As for the \emph{iHS} (see \ref{def.ihs}) and \emph{Rsb}, \emph{xpEHH} is constructed to have an approximately standard Gaussian distribution and may further be transformed into $p_\mathrm{xpEHH}$:
\begin{equation}
p_\mathrm{xpEHH}=-\log_{10}\left(1-2|\Phi\left(\mathrm{xpEHH}\right)-0.5|\right)
\end{equation}
where $\Phi\left(x\right)$ represents the Gaussian cumulative distribution function. Assuming most of the genotyped SNPs behave neutrally (i.e., the genome-wide empirical \emph{xpEHH} distribution is a fair approximation of their corresponding neutral distributions), $p_\mathrm{xpEHH}$ might thus be interpreted as a two-sided P-value (in a $-\log_{10}$ scale) associated to the neutral hypothesis of no selection.
Alternatively, $p_\mathrm{xpEHH}$ might also be computed [@Gautier2011]:
\begin{equation}
p\prime_\mathrm{xpEHH}=-\log_{10}\left(|\Phi\left(\mathrm{xpEHH}\right)|\right)
\end{equation}
$p\prime_\mathrm{xpEHH}$ and $p\prime_\mathrm{xpEHH}$ might then be interpreted as a one-sided P-value (in a $-\log_{10}$ scale) allowing the identification of those sites displaying a significantly high extended haplotype homozygosity in population $pop2$ (represented in the denominator of the corresponding $\mathrm{LRiES}$) relatively to the $pop1$ reference population.

###The function \texttt{ies2xpehh()}
\label{res.xpehh}

The \texttt{ies2xpehh()} function allows to compute \emph{xpEHH} using two data frames containing the \emph{iES} statistics for each of the two populations considered in the same format as the one obtained after running the \texttt{scan\_hh()} function (see \ref{fun.scanhh}).
For instance, to perform a genome scan one might first run for each population \texttt{scan\_hh()} in turn on haplotype data from each chromosome and concatenate the resulting matrices. In the following example, we assume that the haplotype files are named as \texttt{hap\_chr\_i.pop1} and  \texttt{hap\_chr\_i.pop2} where $i$ is the chromosome number (going from 1 to 29), the suffixes pop1 and pop2 indicate the population of origin and the SNP information file is named \texttt{snp.info}. The R code below then generates two data frames (\texttt{wg.res.pop1} and \texttt{wg.res.pop2}) containing the results from all SNPs in the appropriate format to compute \emph{Rsb} with the \texttt{ies2rsb()} function:

```{r,eval=FALSE}
for(i in 1:29){
 hap_file=paste("hap_chr_",i,".pop1",sep="")
 data<-data2haplohh(hap_file="hap_file","snp.info",chr.name=i)
 res<-scan_hh(data)
 if(i==1){wg.res.pop1<-res}else{wg.res.pop1<-rbind(wg.res.pop1,res)}
 hap_file=paste("hap_chr_",i,".pop2",sep="")
 data<-data2haplohh(hap_file="hap_file","snp.info",chr.name=i)
 res<-scan_hh(data)
 if(i==1){wg.res.pop2<-res}else{wg.res.pop2<-rbind(wg.res.pop2,res)}
 }
wg.xpehh<-ies2xpehh(wg.res.pop1,wg.res.pop2)
```
As a matter of illustration, one may consider results from a similar genome scan [@Gautier2011] provided as example data sets and compute for each SNP the \emph{xpEHH} between the CGU and EUT populations as follows:
```{r}
data(wgscan.cgu) ; data(wgscan.eut)
## results from a genome scan (44,057 SNPs) see ?wgscan.eut and ?wgscan.cgu for details
cguVSeut.xpehh<-ies2xpehh(wgscan.cgu,wgscan.eut,"CGU","EUT")
```

The resulting object \texttt{cguVSeut.xpehh} is a data frame with of SNP \emph{xpEHH} (and corresponding P-values assuming \emph{xpEHH} are normally distributed under the neutral hypothesis). Note that either bilateral (default) or unilateral might be performed (\texttt{method} argument). The five first rows of this data frame are displayed below using the following R command:

```{r}
head(cguVSeut.xpehh)
```
###Manhattan plot of the results: the function \texttt{xpehhplot()}

The \texttt{xpehhplot()} function allows to draw a Manhattan plot of the Whole Genome scan results as stored in the data frame produced by the function \texttt{ies2xpehh()}.  Various options are available to modify the graphical display (see \texttt{?xpehhplot}).
As an example, the Figure \ref{fig:xpehhplot} below provides the output of the function \texttt{xpehhplot} for the \emph{xpEHH} computed above across the CGU and EUT populations (see \ref{res.xpehh}).  Figure \ref{fig:xpehhplot} was drawn using the following R code:

```{r,eval=FALSE}
xpehhplot(cguVSeut.xpehh,plot.pval=TRUE)
```

```{r xpehhplot,echo=FALSE,fig.align='center',fig.ext="png",fig.width=16,fig.height=12,fig.lp='fig:',fig.cap='Graphical output for the \\texttt{xpehhplot()} function',fig.pos="!h"}
layout(matrix(1:2,2,1))
plot.pval = TRUE; ylim.scan = 2; pch = 16; cex = 0.75; cex.lab = 1; main = "XPEHH (CGU vs. EUT)"; cex.main = 1.;cex.axis=0.7
data = cguVSeut.xpehh
lst_chrm <- unique(data$CHR) ; nbr_chrm <- length(lst_chrm) ; cum <- rep(0, (nbr_chrm + 1))
if (nbr_chrm > 1) {
  for (i in 1:nbr_chrm) {cum[i + 1] <- tail(data$POSITION[data$CHR == lst_chrm[i]],n = 1)}
}
cum <- cumsum(cum) ; pos <- rep(0, length(data$CHR)) ; pos_labels <- rep(0, nbr_chrm)
for (i in 1:nbr_chrm) {
  pos[data$CHR == lst_chrm[i]] <- data$POSITION[data$CHR ==lst_chrm[i]] + cum[which(lst_chrm == lst_chrm[i])]
  pos_labels[i] <- (cum[i] + cum[i + 1])/2
}
col_chr <- data$CHR
par(mar = c(5, 5, 4, 2) + 0.1)
plot(pos, data[,3], pch = pch, cex = cex, las = 1, col = col_chr,xaxt = "n", xlab = "Chromosome", ylab = expression(italic(XPEHH)), 
     cex.lab = cex.lab, main = main, cex.main = cex.main,cex.axis=cex.axis)
axis(1, at = pos_labels, labels = lst_chrm, las = 1,cex.lab=cex.lab,cex.axis=cex.axis)
abline(h = c(-ylim.scan, ylim.scan), lty = 2)
par(mar = c(5, 5, 4, 2) + 0.1)
plot(pos, data$"-log10(p-value)", pch = pch, cex = cex,las = 1, col = col_chr, xaxt = "n", xlab = "Chromosome", 
         ylab = expression("-" * log[10] * "[1" ~ "-" ~"2" ~ "|" ~ Phi[scriptstyle(italic(XPEHH))] ~"-" ~ 0.5 * "|]"), cex.lab = cex.lab, main = main, cex.main = cex.main,cex.axis=cex.axis)
axis(1, at = pos_labels, labels = lst_chrm, las = 1,cex.lab=cex.lab,cex.axis=cex.axis)
abline(h = c(-ylim.scan, ylim.scan), lty = 2)
```
###\emph{xpEHH} vs. \emph{Rsb} comparison:

A plot of the \emph{xpEHH} against \emph{Rsb} estimates across the CGU and EUT populations (see \ref{res.xpehh} and \ref{res.rsb} respectively) is represented in the Figure \ref{fig:comp} below. This figure was generated using the following R code:

```{r comp,echo=TRUE,fig.align='center',fig.ext="png",out.height="7cm",fig.lp='fig:',fig.cap='Comparison of the XPEHH and Rsb estimates across the CGU and EUT populations',fig.pos="!h"}
plot(cguVSeut.rsb[,3],cguVSeut.xpehh[,3],xlab="Rsb",ylab="XPEHH",pch=16,
    cex=0.5,cex.lab=0.75,cex.axis=0.75)
abline(a=0,b=1,lty=2)
```
##Visual inspection of the standardized scores distribution: the function \texttt{distribplot()}

The \texttt{distribplot} function allows to easily visualize the distributions of the standardized scores (either \emph{iHS}, \emph{Rsb} or \emph{xpEHH}) and compare them to the standard Gaussian distribution. As an example, the Figure \ref{fig:distribplot} below provides the output the function \texttt{distribplot} when considering the \emph{iHS} estimates obtained for the CGU population (see \ref{cgu.ihs}) using the following R code:

```{r,eval=FALSE}
distribplot(ihs.cgu$iHS[,3],xlab="iHS")
```

```{r distribplot,echo=FALSE,fig.align='center',fig.width=7,out.height="7cm",fig.height=7,fig.lp='fig:',fig.cap='Graphical output for the function \\texttt{distribplot}',fig.pos="!h"}
layout(matrix(1:2,2,1))
data=ihs.cgu$iHS[,3]
lty = 1;lwd = 1.5;col = c("blue","red");main = "Genome-wide distribution";xlab = "iHS";cex.main = 1.5;cex.lab = 1.25;qqplot = TRUE
plot(density(data,na.rm = TRUE),main = main,xlab = xlab,col = col[1],lty = lty,lwd = lwd,cex.main = cex.main,cex.lab = cex.lab)
	curve(dnorm,col = col[2],add = TRUE)
	legend("topright",c("Observed","Gaussian"),bty = "n",col = col,lty = lty,lwd = lwd)
	if (qqplot) {
		par(mar = c(5,5,4,2) + 0.1)
		qqnorm(data[!is.na(data)],cex.main = cex.main,cex.lab = cex.lab,pch = 16,cex = 0.75)
     abline(a=0,b=1,lty=2)
    }
```

#Visualizing haplotype structure around a core allele: the function \texttt{bifurcation.diagram()}

The function \texttt{bifurcation.diagram()} function draws haplotype bifurcation diagrams [@Sabeti2002] that allow to better understand the origin of an observed footprints of selection. Such diagrams indeed consist in plotting the breakdown of LD at increasing distances from the core allele at the selected focal SNPs. The root (focal SNP) of each diagram is the core allele and is here identified by a vertical dashed line. 
The diagram is bi-directional, portraying both centromere-proximal and centromere-distal LD. Moving in one direction, each marker is an opportunity for a node; the diagram either divides or not based on whether both or only one allele is present. Thus the breakdown of LD on the core haplotype background is portrayed at progressively longer distances. The thickness of the lines corresponds to the number of samples with the indicated long-distance haplotype. Several options are available to modify the aspect of the plots (see command \texttt{?bifurcation.diagram})
As a matter of illustration, Figure \ref{fig:bifdia} shows the bifurcation diagrams for both the derived and ancestral alleles at the $456^{\text{th}}$ SNP on BTA12 CGU haplotypes. This SNP displayed a strong signal of selection (using both \emph{iHS} and \emph{Rsb} statistics) and is located closed (<5kb) to a strong candidate genes involved in horn development [@Gautier2011]. Figure \ref{fig:bifdia} was obtained with the following R code: 

```{r bifdia,eval=TRUE,fig.align='center',out.height="12cm",fig.width=16,fig.height=12,fig.lp='fig:',fig.cap='Graphical output for the function \\texttt{bifurcation.diagram()}',fig.pos="!h"}
data(haplohh_cgu_bta12)
layout(matrix(1:2,2,1))
bifurcation.diagram(haplohh_cgu_bta12,mrk_foc=456,all_foc=1,nmrk_l=20,nmrk_r=20,
                    main="Bifurcation diagram (RXFP2 SNP on BTA12): Ancestral Allele")
bifurcation.diagram(haplohh_cgu_bta12,mrk_foc=456,all_foc=2,nmrk_l=20,nmrk_r=20,
                    main="Bifurcation diagram (RXFP2 SNP on BTA12): Derived Allele")
```

\clearpage
\newpage

#References