admixture_graph/README.Rmd at master · mailund/admixture_graph · GitHub

1
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
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
---
output:
  md_document:
    variant: markdown_github
---

<!-- README.md is generated from README.Rmd. Please edit that file -->

```{r, echo = FALSE, warning=FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  fig.path = "README-"
)

library(admixturegraph)
suppressPackageStartupMessages(library(dplyr, quietly = TRUE))
suppressPackageStartupMessages(library(ggplot2, quietly = TRUE))
suppressPackageStartupMessages(library(ggthemes, quietly = TRUE))
suppressPackageStartupMessages(library(knitr, quietly = TRUE))
suppressPackageStartupMessages(library(neldermead, quietly = TRUE))
```

# Admixture Graph Manipulation and Fitting

The package provides functionality to analyse and test admixture graphs against the *f* statistics described in the paper [Ancient Admixture in Human History](http://tinyurl.com/o5a4kr4), Patternson *et al.*, Genetics, Vol. 192, 1065--1093, 2012.

The *f* statistics --- *f2*, *f3*, and *f4* --- extract information about correlations between gene frequencies in different populations (or single diploid genome samples), which can be informative about patterns of gene flow between these populations in form of admixture events. If a graph is constructed as a hypothesis for the relationship between the populations, equations for the expected values of the *f* statistics can be extracted, as functions of edge lenghs --- representing genetic drift --- and admixture proportions.

This package provides functions for extracting these equations and for fitting them against computed *f* statistics. It does not currently provide functions for computing the *f* statistics --- for that we refer to the [ADMIXTOOLS](https://github.com/DReichLab/AdmixTools) software package.

## Example

Below is a quick example of how the package can be used. The example uses data
from polar bears and brown bears with a black bear as outgroup and is taken from
[Genomic evidence of geographically widespread effect of gene flow from polar
bears into brown bears](http://onlinelibrary.wiley.com/doi/10.1111/mec.13038/abstract).

The BLK sample is the black bear, the PB sample is a polar bear, and the rest
are brown bears.

I have taken the $f$ statistics from Table 1 in the paper:

```{r}
data(bears)
bears
```

The `D` column is the f4(W,X;Y,Z) statistic and the `Z` column is the $Z$-values
obtained from a blocked jacknife (see Patterson *et al.* for details).

From the statistics we can see that the ABC bears (Adm, Bar and Chi) are closer
related to the polar bears compared to the other brown bears. The paper explains
this with gene flow from polar bears into the ABC bears and going further out
from there, but we can also explain this by several waves of admixture from
ancestral polar bears into brown bears:

```{r graph}
leaves <- c("BLK", "PB",
            "Bar", "Chi1", "Chi2", "Adm1", "Adm2",
            "Denali", "Kenai", "Sweden")
inner_nodes <- c("R", "PBBB",
                 "Adm", "Chi", "BC", "ABC",
                 "x", "y", "z",
                 "pb_a1", "pb_a2", "pb_a3", "pb_a4",
                 "bc_a1", "abc_a2", "x_a3", "y_a4")

edges <- parent_edges(c(edge("BLK", "R"),
                        edge("PB", "pb_a1"),
                        edge("pb_a1", "pb_a2"),
                        edge("pb_a2", "pb_a3"),
                        edge("pb_a3", "pb_a4"),
                        edge("pb_a4", "PBBB"),

                        edge("Chi1", "Chi"),
                        edge("Chi2", "Chi"),
                        edge("Chi", "BC"),
                        edge("Bar", "BC"),
                        edge("BC", "bc_a1"),

                        edge("Adm1", "Adm"),
                        edge("Adm2", "Adm"),

                        admixture_edge("bc_a1", "pb_a1", "ABC", "a"),
                        edge("Adm", "ABC"),

                        edge("ABC", "abc_a2"),
                        admixture_edge("abc_a2", "pb_a2", "x", "b"),

                        edge("Denali", "x"),
                        edge("x", "x_a3"),
                        admixture_edge("x_a3", "pb_a3", "y", "c"),

                        edge("Kenai", "y"),
                        edge("y", "y_a4"),
                        admixture_edge("y_a4", "pb_a4", "z", "d"),

                        edge("Sweden", "z"),

                        edge("z", "PBBB"),
                        edge("PBBB", "R")))

bears_graph <- agraph(leaves, inner_nodes, edges)
plot(bears_graph, show_admixture_labels = TRUE)
```

## Fitting a graph to data

The graph makes predictions on how the *f4* statistics should look. The graph parameters can be fit to observed statistics using the `fit_graph` function:

```{r, warning=FALSE, cache=TRUE, dependson="graph"}
fit <- fit_graph(bears, bears_graph)
fit
```

You can get detailsabout the fit by calling the `summary.agraph_fit` function:

```{r}
summary(fit)
```

You can make a plot of the fit against the data by calling the `plot.agraph_fit` function:

```{r fitted_data}
plot(fit)
```

The plot shows the observed *f4* statistics with error bars (in black) plus the predicted values from the graph.

The result of this is a `ggplot2` object that you can modify by adding `ggplot2` commands in the usual way.

Read the vignette `admixturegraph` for more examples.