---
title: "Evolutionary biology in R (evobiR) tutorial"
author: "Michelle Jonika, Jeffery Demuth, and Richard A. Adams, and Heath Blackmon "
date: "last updated: `r Sys.Date()`"
output:
  rmarkdown::html_vignette:
      fig_caption: yes
      toc: true
vignette: >
  %\VignetteIndexEntry{evobiR tutorial}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

_______
# Installing

The most recent version may be installed from github using the devtools package: 

Installing from github

    library(devtools)
    install_github("coleoguy/evobiR", build_vignettes = TRUE)

```{r, echo=F, tidy=T}
library(evobiR)
```

_______
# Comparative Methods

## AncCond: 
```{}
data(mite.trait)
data(trees.mite)
AncCond(trees[[1]], mite.trait) 
```

This function uses stochastic mapping and ancestral state reconstruction to determine if the derived state of a binary trait originates when a continuous trait has an extreme value.
 
<img src="./images/anccond.png" width="500" height="500"/>

The four steps in the ancestral condition test. A) Ancestral states are estimated for the continuous character assuming a Brownian motion model of evolution. B) Possible evolutionary histories for the discrete trait are generated through stochastic mapping C) Nodes are categorized as either ancestral or derived, and ancestral nodes that subtend an origin of the derived state of the discrete character are annotated. D) Depiction of the null distribution and the observed mean of producing nodes estimated from the data. In this example the producing nodes have a lower continuous value than expected if there is no relationship between the traits.

&nbsp;


 
&nbsp;
 

 
## FuzzyMatch
When assembling data from different sources typos can sometimes cause a loss of perfect matches between trees and datasets. This function helps you find these close matches that can be hand curated to keep as many species as possible in your analysis.

```{r}
data(hym.tree)
names <- c("Pepsis_elegans", "Plagiolepis_alluaudi", "Pheidele_lucreti", "Meliturgula_scriptifronsi", "Andrena_afimbriat")
FuzzyMatch(tree = hym.tree, data = names, max.dist=3)
```

&nbsp;
 

&nbsp;


## SampleTrees

```{r, eval=F}
SampleTrees(trees = system.file("trees.nex", package = "evobiR"), 
            burnin = .1, final.number = 20, format = 'new', prefix = 'sample')
```

This function takes as its input a large collection of trees from a program like MrBayes or Beast and allows the user to select the number of randomly drawn trees they wish to retrieve and to save them in either newick or nexus format.


&nbsp;


&nbsp;


## SuperMatrix

`SuperMatrix(missing = "N", prefix = "DATASET2", save = T)`

This function reads all fasta format alignments in the working directory and constructs a single supermatrix that includes all taxa present in any of the fasta files and inserts missing symbols for taxa that are missing sequences for some loci.  A list with two elements is returned. The first element contains partition data that records the alignment positions of each input fasta file in the combined supermatrix. The second element is a dataframe version of the supermatrix. If the argument save is set to True then both of these files are also saved to the working directory.

&nbsp;

_______
# Population Genetics

## CalcD & WinCalcD
The functions CalcD and CalcPopD are implementations of the algorithm introgression in genomic data. Significance of the D-stat can be calculated either through bootstrapping or jackknifing. Bootstrapping is appropriate for datasets where SNPs are unlinked for instance unmapped RADSeq data. Jackknifing is the appropriate approach when SNPs are potentially in linkage for instance gene alignments or genome alignments. The WinCalcD function is identical to CalcD but performs testing in windows along an alignment.

Durand, Eric Y., et al. Testing for ancient admixture between closely related populations. Molecular biology and evolution 28.8 (2011): 2239-2252.

Eaton, D. A. R., and R. H. Ree. 2013. Inferring phylogeny and introgression using RADseq data: An example from flowering plants (Pedicularis: Orobanchaceae). Syst. Biol. 62:689-706


```{r}
CalcD(alignment = system.file("1.fasta", package = "evobiR"), sig.test = "N")

CalcPopD(alignment = system.file("3.fasta", package = "evobiR"), sig.test = "N")

```


&nbsp;

_______
# Utility/Misc. Functions


## Sliding window

Applies a function within a sliding window of a numeric vector. Both the step size and the window size can be set by the user.  Sliding window analyses are important tools particularly during data exploration.  Often we can find patterns at scales that we might miss otherwise.  As an example lets look at the sunspot data included in R.

```{r, fig.height=5, fig.width=5}
data(sunspot.month)
foo <- as.vector(sunspot.month)
plot(foo, ylab="Number of sunspots", xlab="record number")
```

When we look at the data at this scale the pattern that really catches our attention is the 25 peaks that we see across these 260 years.  This is the well documented 11-year sunspot cycle.  However, our sun has longer cycles that are less evident in this graphing.

```{r, fig.height=5, fig.width=5}
# first lets use the sliding window function to get the number sun spots 

# average over 11 years and do this moving in 1 year steps through time
sunspots <- SlidingWindow(FUN = "mean", 
                          data = foo, 
                          window = 132, 
                          step = 12, 
                          strict = F)

# we repeat this on the years so we can plot agains a sensible x axis
years <- round(SlidingWindow("mean", 
                             data = rep(1749:2013, each=12)[1:3177], 
                             window = 132, 
                             step = 12, 
                             strict=F))
{plot(x=years,y=sunspots, type="l", lwd=3)
abline(v=1810, col="red",lwd=3)
text(y=90, x=1810, "Dalton Min.", col="red",pos=4)}
```

now we can easily spot just how exceptional the Dalton minimum of the early 1800s was.

&nbsp;

## ResSel

This function takes measurements of multiple traits and performs a linear regression and identifies those records with the largest and smallest residual. Originally it was written to perform a regression of horn size on body size allowing for high and low selection lines.  It allows users to choose the trait to select on and the trait to control for.  It also lets the user choose the number of individuals selected (strength of selection).

```{r}
data <- read.csv(file = system.file("horn.beetle.csv", package = "evobiR"))
```
The first column of the data file should contain the identifier i.e. the specimen ID or vial that the measurement is from while the traits should be in the next two columns.

```{r}
data[1:10,]
```

We can then run the residual selection function and it will provide us with both a visual depiction of the data and will return a list with elements (high line and low line providing us with ID numbers of selected individuals.

```{r, fig.height=5, fig.width=5}
ResSel(data = data, traits = c(2,3), percent = 15, identifier = 1, model = "linear")
```


&nbsp;

____
For questions or comments contact [Heath Blackmon](mailto:coleoguy@gmail.com)
