Skip to content

Instantly share code, notes, and snippets.

@dereneaton

dereneaton/R-amova.ipynb

Created September 5, 2017 15:17
Show Gist options
  • Star 0 You must be signed in to star a gist
  • Fork 1 You must be signed in to fork a gist
  • Save dereneaton/a8bafdcb9e1718e4ecd8b2baabbb3d61 to your computer and use it in GitHub Desktop.
Save dereneaton/a8bafdcb9e1718e4ecd8b2baabbb3d61 to your computer and use it in GitHub Desktop.
Download ZIP
R-amova-tutorial
Display the source blob
Display the rendered blob
Raw
R-amova.ipynb
{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# AMOVA in R using the 'poppr' library\n",
"\n",
"Following tutorial from https://grunwaldlab.github.io/Population_Genetics_in_R/AMOVA.html"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Preparing data for poppr analysis explained\n",
"https://grunwaldlab.github.io/Population_Genetics_in_R/Data_Preparation.html"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Installation in R 3.2"
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": [
"## install.packages(\"devtools\")\n",
"## library(devtools)\n",
"## install_github(\"igraph/rigraph\")\n",
"## install.packages(\"phangorn\")\n",
"## install.packages(\"pegas\")\n",
"## intsall.packages(\"poppr\")"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"Loading required package: adegenet\n",
"Loading required package: ade4\n",
"Warning message:\n",
"“replacing previous import by ‘expm::expm’ when loading ‘spdep’”\n",
" /// adegenet 2.0.1 is loaded ////////////\n",
"\n",
" > overview: '?adegenet'\n",
" > tutorials/doc/questions: 'adegenetWeb()' \n",
" > bug reports/feature requests: adegenetIssues()\n",
"\n",
"\n",
"Warning message:\n",
"“replacing previous import by ‘magrittr::%>%’ when loading ‘poppr’”This is poppr version 2.4.1. To get started, type package?poppr\n",
"OMP parallel support: available\n"
]
}
],
"source": [
"## load the poppr library\n",
"library('poppr')"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Load a data set"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"/// GENIND OBJECT /////////\n",
"\n",
" // 187 individuals; 56 loci; 56 alleles; size: 66.5 Kb\n",
"\n",
" // Basic content\n",
" @tab: 187 x 56 matrix of allele counts\n",
" @loc.n.all: number of alleles per locus (range: 56-56)\n",
" @ploidy: ploidy of each individual (range: 2-2)\n",
" @type: PA\n",
" @call: old2new_genind(object = x, donor = new(class(x)))\n",
"\n",
" // Optional content\n",
" @pop: population of each individual (group size range: 90-97)\n",
" @other: a list containing: population_hierarchy \n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## load the example data set\n",
"data(\"Aeut\")\n",
"\n",
"## this is a summary of the Aeut data object\n",
"Aeut"
]
},
{
"cell_type": "code",
"execution_count": 41,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<thead><tr><th></th><th scope=col>L01</th><th scope=col>L02</th><th scope=col>L03</th><th scope=col>L04</th><th scope=col>L05</th><th scope=col>L06</th><th scope=col>L07</th><th scope=col>L08</th><th scope=col>L09</th><th scope=col>L10</th><th scope=col>⋯</th><th scope=col>L47</th><th scope=col>L48</th><th scope=col>L49</th><th scope=col>L50</th><th scope=col>L51</th><th scope=col>L52</th><th scope=col>L53</th><th scope=col>L54</th><th scope=col>L55</th><th scope=col>L56</th></tr></thead>\n",
"<tbody>\n",
"\t<tr><th scope=row>001</th><td>1</td><td>0</td><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td><td>0</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>1</td></tr>\n",
"\t<tr><th scope=row>002</th><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td></tr>\n",
"\t<tr><th scope=row>003</th><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
"\t<tr><th scope=row>004</th><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
"\t<tr><th scope=row>005</th><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td></tr>\n",
"\t<tr><th scope=row>006</th><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>1</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>0</td><td>1</td><td>1</td></tr>\n",
"\t<tr><th scope=row>007</th><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
"\t<tr><th scope=row>008</th><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>0</td><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
"\t<tr><th scope=row>009</th><td>1</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td></tr>\n",
"\t<tr><th scope=row>010</th><td>1</td><td>1</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>0</td><td>1</td><td>0</td><td>⋯</td><td>0</td><td>0</td><td>1</td><td>1</td><td>1</td><td>1</td><td>1</td><td>0</td><td>1</td><td>0</td></tr>\n",
"</tbody>\n",
"</table>\n"
],
"text/latex": [
"\\begin{tabular}{r|llllllllllllllllllllllllllllllllllllllllllllllllllllllll}\n",
" & L01 & L02 & L03 & L04 & L05 & L06 & L07 & L08 & L09 & L10 & ⋯ & L47 & L48 & L49 & L50 & L51 & L52 & L53 & L54 & L55 & L56\\\\\n",
"\\hline\n",
"\t001 & 1 & 0 & 1 & 0 & 1 & 0 & 0 & 1 & 1 & 0 & ⋯ & 0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 & 0 & 1\\\\\n",
"\t002 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 1 & 1 & 1 & ⋯ & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 0 & 0 & 1\\\\\n",
"\t003 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & ⋯ & 0 & 0 & 1 & 0 & 1 & 1 & 0 & 0 & 0 & 0\\\\\n",
"\t004 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & ⋯ & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\\\\n",
"\t005 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & ⋯ & 0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 & 1 & 0\\\\\n",
"\t006 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 1 & ⋯ & 0 & 0 & 1 & 1 & 1 & 1 & 0 & 0 & 1 & 1\\\\\n",
"\t007 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & ⋯ & 0 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 0\\\\\n",
"\t008 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & ⋯ & 0 & 0 & 1 & 0 & 1 & 1 & 0 & 0 & 0 & 0\\\\\n",
"\t009 & 1 & 0 & 1 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & ⋯ & 0 & 0 & 1 & 0 & 0 & 1 & 0 & 0 & 0 & 0\\\\\n",
"\t010 & 1 & 1 & 0 & 0 & 0 & 0 & 0 & 0 & 1 & 0 & ⋯ & 0 & 0 & 1 & 1 & 1 & 1 & 1 & 0 & 1 & 0\\\\\n",
"\\end{tabular}\n"
],
"text/markdown": [
"\n",
"| <!--/--> | L01 | L02 | L03 | L04 | L05 | L06 | L07 | L08 | L09 | L10 | ⋯ | L47 | L48 | L49 | L50 | L51 | L52 | L53 | L54 | L55 | L56 | \n",
"|---|---|---|---|---|---|---|---|---|---|\n",
"| 001 | 1 | 0 | 1 | 0 | 1 | 0 | 0 | 1 | 1 | 0 | ⋯ | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 0 | 1 | \n",
"| 002 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 1 | 1 | 1 | ⋯ | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | \n",
"| 003 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | ⋯ | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | \n",
"| 004 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | ⋯ | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | \n",
"| 005 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | ⋯ | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 0 | \n",
"| 006 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | ⋯ | 0 | 0 | 1 | 1 | 1 | 1 | 0 | 0 | 1 | 1 | \n",
"| 007 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | ⋯ | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0 | \n",
"| 008 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | ⋯ | 0 | 0 | 1 | 0 | 1 | 1 | 0 | 0 | 0 | 0 | \n",
"| 009 | 1 | 0 | 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | ⋯ | 0 | 0 | 1 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | \n",
"| 010 | 1 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | ⋯ | 0 | 0 | 1 | 1 | 1 | 1 | 1 | 0 | 1 | 0 | \n",
"\n",
"\n"
],
"text/plain": [
" L01 L02 L03 L04 L05 L06 L07 L08 L09 L10 ⋯ L47 L48 L49 L50 L51 L52 L53 L54\n",
"001 1 0 1 0 1 0 0 1 1 0 ⋯ 0 0 1 1 1 1 0 0 \n",
"002 1 0 1 0 0 0 0 1 1 1 ⋯ 0 0 1 1 1 1 1 0 \n",
"003 1 0 1 0 0 0 0 0 1 0 ⋯ 0 0 1 0 1 1 0 0 \n",
"004 1 0 1 0 0 0 0 0 1 0 ⋯ 0 0 1 0 0 1 0 0 \n",
"005 1 0 1 0 0 0 0 0 1 0 ⋯ 0 0 1 1 1 1 0 0 \n",
"006 1 0 1 0 0 0 0 0 1 1 ⋯ 0 0 1 1 1 1 0 0 \n",
"007 1 0 1 0 0 0 0 0 1 0 ⋯ 0 0 1 0 0 0 0 0 \n",
"008 1 0 1 0 0 0 0 0 1 0 ⋯ 0 0 1 0 1 1 0 0 \n",
"009 1 0 1 0 0 0 0 0 1 0 ⋯ 0 0 1 0 0 1 0 0 \n",
"010 1 1 0 0 0 0 0 0 1 0 ⋯ 0 0 1 1 1 1 1 0 \n",
" L55 L56\n",
"001 0 1 \n",
"002 0 1 \n",
"003 0 0 \n",
"004 0 0 \n",
"005 1 0 \n",
"006 1 1 \n",
"007 0 0 \n",
"008 0 0 \n",
"009 0 0 \n",
"010 1 0 "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## this is the genotype data (binary in this case)\n",
"head(Aeut$tab, 10)"
]
},
{
"cell_type": "code",
"execution_count": 4,
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
"<table>\n",
"<thead><tr><th scope=col>Pop_Subpop</th><th scope=col>Pop</th><th scope=col>Subpop</th></tr></thead>\n",
"<tbody>\n",
"\t<tr><td>Athena_1</td><td>Athena </td><td>1 </td></tr>\n",
"\t<tr><td>Athena_1</td><td>Athena </td><td>1 </td></tr>\n",
"\t<tr><td>Athena_1</td><td>Athena </td><td>1 </td></tr>\n",
"\t<tr><td>Athena_1</td><td>Athena </td><td>1 </td></tr>\n",
"\t<tr><td>Athena_1</td><td>Athena </td><td>1 </td></tr>\n",
"\t<tr><td>Athena_1</td><td>Athena </td><td>1 </td></tr>\n",
"</tbody>\n",
"</table>\n"
],
"text/latex": [
"\\begin{tabular}{r|lll}\n",
" Pop\\_Subpop & Pop & Subpop\\\\\n",
"\\hline\n",
"\t Athena\\_1 & Athena & 1 \\\\\n",
"\t Athena\\_1 & Athena & 1 \\\\\n",
"\t Athena\\_1 & Athena & 1 \\\\\n",
"\t Athena\\_1 & Athena & 1 \\\\\n",
"\t Athena\\_1 & Athena & 1 \\\\\n",
"\t Athena\\_1 & Athena & 1 \\\\\n",
"\\end{tabular}\n"
],
"text/markdown": [
"\n",
"Pop_Subpop | Pop | Subpop | \n",
"|---|---|---|---|---|---|\n",
"| Athena_1 | Athena | 1 | \n",
"| Athena_1 | Athena | 1 | \n",
"| Athena_1 | Athena | 1 | \n",
"| Athena_1 | Athena | 1 | \n",
"| Athena_1 | Athena | 1 | \n",
"| Athena_1 | Athena | 1 | \n",
"\n",
"\n"
],
"text/plain": [
" Pop_Subpop Pop Subpop\n",
"1 Athena_1 Athena 1 \n",
"2 Athena_1 Athena 1 \n",
"3 Athena_1 Athena 1 \n",
"4 Athena_1 Athena 1 \n",
"5 Athena_1 Athena 1 \n",
"6 Athena_1 Athena 1 "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## format the 'other$pop_hier' data into 'strata' and show it\n",
"strata(Aeut) <- data.frame(other(Aeut)$population_hierarchy)\n",
"head(Aeut$strata)"
]
},
{
"cell_type": "code",
"execution_count": 5,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\n",
"This is a genclone object\n",
"-------------------------\n",
"Genotype information:\n",
"\n",
" 119 original multilocus genotypes \n",
" 187 diploid individuals\n",
" 56 dominant loci\n",
"\n",
"Population information:\n",
"\n",
" 3 strata - Pop_Subpop, Pop, Subpop\n",
" 2 populations defined - Athena, Mt. Vernon"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## make in genclone object\n",
"Aeut <- as.genclone(Aeut)\n",
"Aeut"
]
},
{
"cell_type": "code",
"execution_count": 6,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"\n",
" Athena Mt. Vernon \n",
" 97 90 "
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## how are individuals distributed among the populations?\n",
"table(strata(Aeut, ~Pop))"
]
},
{
"cell_type": "code",
"execution_count": 7,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
" Subpop\n",
"Pop 1 2 3 4 5 6 7 8 9 10\n",
" Athena 9 12 10 13 10 5 11 8 10 9\n",
" Mt. Vernon 10 6 8 12 17 12 12 13 0 0"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"## how are they distributed among the subpopulations\n",
"table(strata(Aeut, ~Pop/Subpop, combine = FALSE)) "
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Run the AMOVA"
]
},
{
"cell_type": "code",
"execution_count": 9,
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"\n",
" No missing values detected.\n",
"\n",
"\n",
" No missing values detected.\n",
"\n"
]
}
],
"source": [
"Aeutamova <- poppr.amova(Aeut, ~Pop/Subpop)\n",
"Aeutamovacc <- poppr.amova(Aeut, ~Pop/Subpop, clonecorrect = TRUE)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"$call\n",
"ade4::amova(samples = xtab, distances = xdist, structures = xstruct)\n",
"\n",
"$results\n",
" Df Sum Sq Mean Sq\n",
"Between Pop 1 1051.2345 1051.234516\n",
"Between samples Within Pop 16 273.4575 17.091091\n",
"Within samples 169 576.5059 3.411277\n",
"Total 186 1901.1979 10.221494\n",
"\n",
"$componentsofcovariance\n",
" Sigma %\n",
"Variations Between Pop 11.063446 70.006786\n",
"Variations Between samples Within Pop 1.328667 8.407483\n",
"Variations Within samples 3.411277 21.585732\n",
"Total variations 15.803391 100.000000\n",
"\n",
"$statphi\n",
" Phi\n",
"Phi-samples-total 0.7841427\n",
"Phi-samples-Pop 0.2803128\n",
"Phi-Pop-total 0.7000679\n"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"Aeutamova"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"collapsed": true
},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "R 3.2",
"language": "R",
"name": "ir32"
},
"language_info": {
"codemirror_mode": "r",
"file_extension": ".r",
"mimetype": "text/x-r-source",
"name": "R",
"pygments_lexer": "r",
"version": "3.2.3"
}
},
"nbformat": 4,
"nbformat_minor": 2
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment