An introduction to adegenet 2.1.

6
*
Thibaut Jombart
Imperial College London
MRC Centre for Outbreak Analysis and Modelling

June 5, 2022

Abstract
This vignette provides an introductory tutorial to the adegenet package [3] for
the R software [11]. This package implements tools to handle, analyse and simulate
genetic data. Originally developped for multiallelic, codominant markers such as
microsatellites, adegenet now also handles dominant markers, allows for any ploidy in
the data, handles SNPs and sequence data, and implements a memory-efficient storage
for genome-wide SNP data. This tutorial provides an overview of adegenet’s basic
functionalities. Since adegenet 1.4-0, this tutorial is no longer distributed as a package
vignette. Also note that adegenet has undergone substantial changes with version
2.0.0, including a reform of the data structure and new accessors, all documented in
this tutorial.

* tjombart@[Link]

1
Contents
1 Introduction 3

2 Getting started 4
2.1 Installing the package - stable version . . . . . . . . . . . . . . . . . . . . . . 4
2.2 Installing the package - devel version . . . . . . . . . . . . . . . . . . . . . . 5
2.3 Getting help in R . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.4 Asking help on a forum . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6
2.5 Bug report, feature requests, contributions: we are all one! . . . . . . . . . . 7

3 Object classes 8
3.1 genind objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
3.2 genpop objects . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 10
3.3 Using accessors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11

4 Importing/exporting data 15
4.1 Importing data from GENETIX, STRUCTURE, FSTAT, Genepop . . . . . 15
4.2 Importing data from other software . . . . . . . . . . . . . . . . . . . . . . . 15
4.3 Handling presence/absence data . . . . . . . . . . . . . . . . . . . . . . . . . 17
4.4 SNPs data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4.5 Extracting polymorphism from DNA sequences . . . . . . . . . . . . . . . . 22
4.6 Extracting polymorphism from proteic sequences . . . . . . . . . . . . . . . 26
4.7 Using genind/genpop constructors . . . . . . . . . . . . . . . . . . . . . . . 30
4.8 Exporting data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32

5 Basics of data analysis 33

5.1 Manipulating the data . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33
5.2 Using summaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
5.3 Testing for Hardy-Weinberg equilibrium . . . . . . . . . . . . . . . . . . . . 42
5.4 Measuring and testing population structure (a.k.a F statistics) . . . . . . . . 42
5.5 Estimating inbreeding . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44

6 Multivariate analysis 49
6.1 General overview . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49
6.2 Performing a Principal Component Analysis on genind objects . . . . . . . . 51
6.3 Performing a Correspondance Analysis on genpop objects . . . . . . . . . . . 58

7 Spatial analysis 63
7.1 Isolation by distance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63
7.1.1 Testing isolation by distance . . . . . . . . . . . . . . . . . . . . . . . 63
7.1.2 Cline or distant patches? . . . . . . . . . . . . . . . . . . . . . . . . . 65
7.2 Using Monmonier’s algorithm to define genetic boundaries . . . . . . . . . . 67

8 Simulating hybridization 77

2
1 Introduction
This tutorial introduces some basic functionalities of the adegenet package for R [11]. The
purpose of this package is to provide tools for handling, analysing and simulating genetic
data, with an emphasis on multivariate approaches and exploratory methods. Standard
multivariate analyses are implemented in the ade4 package [2], of which adegenet was
originally an extension. However, the package has since grown methods of its own such
as the Discriminant Analysis of Principal Components (DAPC, [7]), the spatial Principal
Components Analysis (sPCA, [4]), or the SeqTrack algorithm [5]. In this tutorial, we
introduce the main data structures, show how to import data into adegenet, and cover some
basic population genetics and multivariate analysis.

Other tutorials are available via the command adegenetTutorial:

ˆ adegenetTutorial("spca"): tutorial on the sPCA

ˆ adegenetTutorial("dapc"): tutorial on the DAPC

ˆ adegenetTutorial("genomics"): tutorial on handling large SNP datasets using

genlight objects

ˆ adegenetTutorial("strata"): tutorial on handling stratified population data

3
2 Getting started
2.1 Installing the package - stable version
Before going further, we shall make sure that adegenet is well installed on the computer. The
current version of the package is 2.1.6. Make sure you have a recent version of R (≥ 3.2.1)
by typing:

[Link]

## [1] "R version 4.2.0 (2022-04-22)"

Then, install adegenet with dependencies using:

[Link]("adegenet", dep=TRUE)

We can now load the package alongside other useful packages:

library("ape")
library("pegas")
library("seqinr")
library("ggplot2")

library("adegenet")

## Loading required package: ade4

##
## Attaching package: ’ade4’
## The following object is masked from ’package:pegas’:
##
## amova
##
## /// adegenet 2.1.6 is loaded ////////////
##
## > overview: ’?adegenet’
## > tutorials/doc/questions: ’adegenetWeb()’
## > bug reports/feature requests: adegenetIssues()

If at some point you are unsure about the version of the package, you can check it using:

packageDescription("adegenet", fields = "Version")

## [1] "2.1.6"

adegenet version should read 2.1.6.

4
2.2 Installing the package - devel version
The development of adegenet is hosted on github:
[Link]
You can install this version using the package devtools and the following commands:

library("devtools")
install_github("thibautjombart/adegenet")
library("adegenet")

The development version may implement new features and fix known issues. However, it
may also occasionally be broken, as this is our working copy of the project. Usual disclaimers
apply here: this package is provided with no warranty, etc. If unsure, use the stable version.

2.3 Getting help in R

There are several ways of getting information about R in general, and about adegenet in
particular. The function [Link] is used to look for help on a given topic. For instance:

[Link]("Hardy-Weinberg")

replies that there is a function [Link] in the adegenet package, and other
similar functions in genetics and pegas. To get help for a given function, use ?foo where foo
is the function of interest. For instance (quotes and parentheses can be removed):

?spca

will open up the manpage of the spatial principal component analysis [4]. At the end
of a manpage, an ‘example’ section often shows how to use a function. This can be copied
and pasted to the console, or directly executed from the console using example. For further
questions concerning R, the function RSiteSearch is a powerful tool for making online
researches using keywords in R’s archives (mailing lists and manpages).

adegenet has a few extra documentation sources. Information can be found from
the website ([Link] in the ‘documents’ section,
including several tutorials and a manual which compiles all manpages of the package, and a
dedicated mailing list with searchable archives. To open the website from R, use:

adegenetWeb()

The same can be done for tutorials, using adegenetTutorial (see manpage to choose the
tutorial to open). You will also find an overview of the main functionalities of the package
typing:

5
?adegenet

Note that you can also browse help pages as html pages, using:

[Link]()

To go to the adegenet page, click ‘packages’, ‘adegenet’, and ‘adegenet-package’.

2.4 Asking help on a forum

Several mailing lists are available to find different kinds of information on R. adegenet has
its own dedicated forum/mailing list: adegenet-forum@[Link]
To avoid spam, this list is filtered; subscription is recommended, and can be done at:
[Link]

Posting questions on R forums can sometimes be a traumatic experience, and we are

trying to avoid this as much as possible on the adegenet forum. To this end, the following
points are worth keeping in mind:

ˆ read the doc first: manpages and tutorials take an awful long time to write and
maintain; make sure your answer is not in an obvious place before asking a question;
pretending to have read all the available doc while you have not even looked at the
basics tutorial is a clever, yet often unsuccessful strategy.

ˆ search the archives: adegenet forum has searchable archives (see the adegenet
website); your answer may be there already, so it is worth checking.

ˆ give us info: you tried something, it is not working.. give us some information: what
version of adegenet are you using, what commands did you enter and what was the
output, etc.

ˆ avoid personal messages: the adegenet forum has plenty of advantages: several
people are likely to reply and participate in the conversation, answers are generally
faster, and all of this is archived and searchable. Please do not email the developers
directly, unless you need to discuss confidential matters.

ˆ short answers are okay: some answers will be short. Do not take them as rude, or
think people are upset: answering questions on a forum is a time-consuming activity
and the reward for it is low. Sometimes the best answer will be pointing to relevant
documentation, e.g. “Please check ?xvalDapc”. If you get this, we (most likely) still
like you.

The adegenet forum is not the only forum that might be relevant. Others include:

ˆ R-sig-genetics: genetics in R.
[Link]

6
 ˆ R-sig-phylo: phylogenetics in R.
[Link]

ˆ R-help: general questions about R.

[Link]

ˆ [Link]: general questions about R.

[Link]

Please avoid double-posting (it is often considered to be rude).

2.5 Bug report, feature requests, contributions: we are all one!

Free software evolves through interactions between communities of developers and users.
This is especially true in R, where these two communities are very much intermingled. In
short: we are all contributors! These contributions include:

ˆ asking a question: see section above

ˆ asking for a new feature: something useful is missing, and you think it will be useful
to others? Say it! Post a feature request using github’s issues:
[Link]

ˆ reporting a possible bug: bugs are rare, but if you think you found one, post it as
an issue on github:
[Link]

ˆ contributions: github makes contributions very easy; fork the project, make the
changes you want, and when you are happy and the package passes the checks, send a
pull request; for more on this go to the github page:
[Link] And please, remember to add
yourself as a contributor in the DESCRIPTION and relevant manpages!

7
3 Object classes
Two main classes of objects are used for storing genetic marker data, depending on the level
at which the genetic information is considered: genind is used for individual genotypes,
whereas genpop is used for alleles numbers counted by populations. Note that the term
’population’, here and later, is employed in a broad sense: it simply refers to any grouping
of individuals. The specific class genlight is used for storing large genome-wide SNPs data.
See the genomics tutorial for more information on this topic.

3.1 genind objects

These objects store genetic data at an individual level, plus various meta-data (e.g. group
membership). genind objects can be obtained by reading data files from other software,
from a [Link] of genotypes, by conversion from a table of allele counts, or even from
aligned DNA or proteic sequences (see ’importing data’). Here, we introduce this class using
the dataset nancycats, which is already stored as a genind object:

data(nancycats)
[Link](nancycats)

## [1] TRUE

nancycats

## /// GENIND OBJECT /////////

##
## // 237 individuals; 9 loci; 108 alleles; size: 150.5 Kb
##
## // Basic content
## @tab: 237 x 108 matrix of allele counts
## @[Link]: number of alleles per locus (range: 8-18)
## @[Link]: locus factor for the 108 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: genind(tab = truenames(nancycats)$tab, pop = truenames(nancycats)$pop)
##
## // Optional content
## @pop: population of each individual (group size range: 9-23)
## @other: a list containing: xy

A genind object is formal S4 object with several slots, accessed using the ’@’ operator
(see class?genind). Note that the ’$’ is also implemented for adegenet objects, so that
slots can be accessed as if they were components of a list.

genind objects possess the following slots:

8
 ˆ tab: a matrix of alleles counts (individuals in rows, alleles in columns).

ˆ [Link]: the number of alleles in each marker.

ˆ [Link]: a factor indicating which columns in @tab correspond to which marker.

ˆ [Link]: a list of allele names for each locus.

ˆ ploidy: a vector of integer indicating the ploidy of each individual; constant ploidy is
assumed across different loci of a single individual, but different individuals can have
different ploidy.

ˆ type: a character string indicating whether the marker is codominant (codom) or

presence/absence (PA).

ˆ call: the matched call, i.e. command used to create the object.

ˆ pop: (optional) a factor storing group membership of the individuals; when present,
method needing a population information will automatically use this if nothing else is
provided.

ˆ strata: (optional) a [Link] storing hierarchical structure of individuals (see

dedicated tutorial).

ˆ hierarchy: (optional) a formula defining the hierarchical structure of individuals (see

dedicated tutorial).

ˆ other: (optional) a list storing optional information.

Slots can be accessed using ’@’ or ’$’, although it is recommended to use accessors to
retrieve and change slot values (see section ’Using accessors’).

The main slot of a genind is the table of allelic counts @tab, with individuals in rows and
alleles in columns. For instance:

nancycats[10:18, loc="fca8"]@tab

## fca8.117 fca8.119 fca8.121 fca8.123 fca8.127 fca8.129 fca8.131 fca8.133 fca8.135

## N224 0 0 0 0 0 0 0 0 2
## N7 0 0 0 0 0 0 0 0 0
## N141 0 0 0 0 0 1 0 1 0
## N142 0 0 0 0 0 1 0 1 0
## N143 0 0 0 0 0 0 0 2 0
## N144 0 0 0 0 0 0 1 0 1
## N145 0 0 0 0 0 1 0 0 1
## N146 0 0 0 0 0 1 0 1 0
## N147 0 0 0 0 0 1 0 0 1

9
 Individual ’N224’ is an homozygote for the allele 135 at the locus ’fca8’, while N141” is
an heterozygote with alleles 129/133. Note that as of adegenet 2.0.0, this table is no storing
data as relative frequencies. These can still be obtained using the accessor ’tab’. The
particular case of presence/absence data is described in a dedicated section (see ’Handling
presence/absence data’). As of version 2.0.0 of adegenet, the slots @strata and @hierarchy
implement hierarchical population structures. See dedicated tutorial for more on this topic.

Note that objects can be regenerated using the matched call stored in genind objects, i.e.
the instruction that created it. For instance:

obj <- [Link]([Link]("files/[Link]",package="adegenet"))

##
## Converting data from GENETIX to a genind object...
##
## ...done.

obj$call

## [Link](file = [Link]("files/[Link]", package = "adegenet"))

toto <- eval(obj$call)

##
## Converting data from GENETIX to a genind object...
##
## ...done.

identical(obj,toto)

## [1] TRUE

3.2 genpop objects

These objects store genetic data at a population level, plus various meta-data. Their struture
is nearly identical to genind objects, only simpler as they no longer store group membership
for individuals. genpop objects are created from genind objects using genind2genpop:

data(nancycats)
catpop <- genind2genpop(nancycats)

##
## Converting data from a genind to a genpop object...
##
## ...done.

10
catpop

## /// GENPOP OBJECT /////////

##
## // 17 populations; 9 loci; 108 alleles; size: 31 Kb
##
## // Basic content
## @tab: 17 x 108 matrix of allele counts
## @[Link]: number of alleles per locus (range: 8-18)
## @[Link]: locus factor for the 108 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: genind2genpop(x = nancycats)
##
## // Optional content
## @other: a list containing: xy

As in genind objects, data are stored as numbers of alleles, but this time for populations
(here, cat colonies):

tab(catpop)[1:5, 1:10] # using the accessor for the allele table

## fca8.117 fca8.119 fca8.121 fca8.123 fca8.127 fca8.129 fca8.131 fca8.133 fca8.135

## P01 0 0 0 0 0 0 0 2 9
## P02 0 0 0 0 0 10 9 8 14
## P03 0 0 0 4 0 0 0 0 1
## P04 0 0 0 3 0 0 0 1 7
## P05 0 0 0 1 0 0 0 0 7

3.3 Using accessors

One advantage of formal (S4) classes is that they allow for interacting simply with possibly
complex objects. This is made possible by using accessors, i.e. functions that extract
information from an object, rather than accessing the slots directly. Another advantage of
this approach is that as long as accessors remain identical on the user’s side, the internal
structure of an object may change from one release to another without generating errors in
old scripts. Although genind and genpop objects are fairly simple, we recommend using
accessors whenever possible to access or modify their content.

Available accessors are:

ˆ nInd: returns the number of individuals in the object; only for genind.

ˆ nLoc: returns the number of loci.

11
 ˆ nAll: returns the number of alleles for each locus.

ˆ nPop: returns the number of populations.

ˆ tab: returns a table of allele numbers, or frequencies (if requested), with optional
replacement of missing values; replaces the former accessor ’truenames’.

ˆ indNames† : returns/sets labels for individuals; only for genind.

ˆ locNames† : returns/sets labels for loci.

ˆ alleles† : returns/sets alleles.

ˆ ploidy† : returns/sets ploidy of the individuals; when setting values, a single value can
be provided, in which case constant ploidy is assumed.

ˆ pop† : returns/sets a factor grouping individuals; only for genind.

ˆ strata† : returns/sets data defining strata of individuals; only for genind.

ˆ hier† : returns/sets hierarchical groups of individuals; only for genind.

ˆ other† : returns/sets misc information stored as a list.

†
where indicates that a replacement method is available using <-; for instance:

head(indNames(nancycats),10)

## [1] "N215" "N216" "N217" "N218" "N219" "N220" "N221" "N222" "N223" "N224"

indNames(nancycats) <- paste("cat", 1:nInd(nancycats),sep=".")

head(indNames(nancycats),10)

## [1] "cat.1" "cat.2" "cat.3" "cat.4" "cat.5" "cat.6" "cat.7" "cat.8" "cat.9"

Some accessors such as locNames may have specific options; for instance:

locNames(nancycats)

## [1] "fca8" "fca23" "fca43" "fca45" "fca77" "fca78" "fca90" "fca96" "fca37"

returns the names of the loci, while:

temp <- locNames(nancycats, withAlleles=TRUE)

head(temp, 10)

## [1] "fca8.117" "fca8.119" "fca8.121" "fca8.123" "fca8.127" "fca8.129" "fca8.131" "fc

12
 returns the names of the alleles in the form ’[Link]’.

The slot ’pop’ can be retrieved and set using pop:

obj <- nancycats[sample(1:50,10)]

pop(obj)

## [1] P04 P02 P02 P04 P01 P03 P02 P02 P02 P02
## Levels: P01 P02 P03 P04

pop(obj) <- rep("newPop",10)

pop(obj)

## [1] newPop newPop newPop newPop newPop newPop newPop newPop newPop newPop
## Levels: newPop

Accessors make things easier. For instance, when setting new names for loci, the columns
of @tab are renamed automatically:

head(colnames(tab(obj)),20)

## [1] "fca8.117" "fca8.119" "fca8.121" "fca8.123" "fca8.127" "fca8.129" "fca8.13

## [15] "fca8.147" "fca8.149" "fca23.128" "fca23.130" "fca23.132" "fca23.136"

locNames(obj)

## [1] "fca8" "fca23" "fca43" "fca45" "fca77" "fca78" "fca90" "fca96" "fca37"

locNames(obj)[1] <- "newLocusName"

locNames(obj)

## [1] "newLocusName" "fca23" "fca43" "fca45" "fca77" "fca78

head(colnames(tab(obj)),20)

## [1] "newLocusName.117" "newLocusName.119" "newLocusName.121" "newLocusName.123" "new

## [10] "newLocusName.137" "newLocusName.139" "newLocusName.141" "newLocusName.143" "new
## [19] "fca23.132" "fca23.136"

An additional advantage of using accessors is they are most of the time safer to use. For
instance, pop<- will check the length of the new group membership vector against the data,
and complain if there is a mismatch. It also converts the provided replacement to a factor,
while the command:

13
obj@pop <- rep("newPop",10)

## Error in (function (cl, name, valueClass) : assignment of an object of

class "character" is not valid for @’pop’ in an object of class "genind";
is(value, "factorOrNULL") is not TRUE

generates an error (since replacement is not a factor).

14
4 Importing/exporting data
4.1 Importing data from GENETIX, STRUCTURE, FSTAT,
Genepop
Data can be read from the software GENETIX (extension .gtx), STRUCTURE (.str or
.stru), FSTAT (.dat) and Genepop (.gen) files, using the corresponding read function:
[Link], [Link], [Link], and [Link]. These functions take
as main argument the path (as a string of characters) to an input file, and produce a genind
object. Alternatively, one can use the function import2genind which detects a file format
from its extension and uses the appropriate routine. For instance:

obj1 <- [Link]([Link]("files/[Link]",package="adegenet"))

##
## Converting data from GENETIX to a genind object...
##
## ...done.

obj2 <- import2genind([Link]("files/[Link]", package="adegenet"))

##
## Converting data from GENETIX to a genind object...
##
## ...done.

[Link](obj1,obj2)

## [1] "Attributes: < Component \"call\": target, current do not match when deparsed >"

The only difference between obj1 and obj2 is their call (which is normal as they were obtained
from different command lines).

4.2 Importing data from other software

Raw genetic markers data are often stored as tables with individuals in row and markers
in column, where each entry is a character string coding the alleles possessed at one locus.
Such data are easily imported into R as a [Link], using for instance [Link] for
text files or [Link] for comma-separated text files. Then, the obtained [Link] can
be converted into a genind object using df2genind.

There are only a few pre-requisite the data should meet for this conversion to be
possible. The easiest and clearest way of coding data is using a separator between alleles.
For instance, ”80/78”, ”80—78”, or ”80,78” are different ways of coding a genotype at a
microsatellite locus with alleles ’80’ and 78”. Note that for haploid data, no separator shall
be used. The only contraint when using a separator is that the same separator is used in all

15
the dataset. There are no contraints as to i) the type of separator used or ii) the ploidy of
the data. These parameters can be set in df2genind through arguments sep and ploidy,
respectively.

Alternatively, no separator may be used provided a fixed number of characters is used

to code each allele. For instance, in a diploid organism, ”0101” is an homozygote 1/1 while
”1209” is a heterozygote 12/09 in a two-character per allele coding scheme. In a tetraploid
system with one character per allele, ”1209” will be understood as 1/2/0/9.
Here, we provide an example using randomly generated tetraploid data and no separator.

temp <- lapply(1:30, function(i) sample(1:9, 4, replace=TRUE))

temp <- sapply(temp, paste, collapse="")
temp <- matrix(temp, nrow=10, dimnames=list(paste("ind",1:10), paste("loc",1:3)))
temp

## loc 1 loc 2 loc 3

## ind 1 "9833" "8832" "1489"
## ind 2 "1795" "6576" "9156"
## ind 3 "1621" "7778" "5266"
## ind 4 "2291" "4868" "2345"
## ind 5 "7165" "4794" "1243"
## ind 6 "8545" "3418" "6814"
## ind 7 "9973" "8888" "5927"
## ind 8 "4492" "6989" "7585"
## ind 9 "6218" "5528" "3668"
## ind 10 "1664" "4542" "1293"

obj <- df2genind(temp, ploidy=4, sep="")

obj

## /// GENIND OBJECT /////////

##
## // 10 individuals; 3 loci; 27 alleles; size: 10.2 Kb
##
## // Basic content
## @tab: 10 x 27 matrix of allele counts
## @[Link]: number of alleles per locus (range: 9-9)
## @[Link]: locus factor for the 27 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 4-4)
## @type: codom
## @call: df2genind(X = temp, sep = "", ploidy = 4)
##
## // Optional content
## - empty -

16
obj is a genind containing the same information, but recoded as a matrix of allele numbers
($tab slot). We can check that the conversion was exact by converting back the object into
a table of character strings (function genind2df):

genind2df(obj, sep="|")

## loc 1 loc 2 loc 3

## ind 1 9|8|3|3 8|8|3|2 1|4|8|9
## ind 2 9|1|7|5 6|6|5|7 1|9|5|6
## ind 3 1|1|6|2 8|7|7|7 5|6|6|2
## ind 4 9|1|2|2 8|8|6|4 4|5|2|3
## ind 5 1|7|5|6 7|4|4|9 1|4|2|3
## ind 6 8|5|5|4 8|3|4|1 1|4|8|6
## ind 7 9|9|3|7 8|8|8|8 9|5|2|7
## ind 8 9|2|4|4 8|6|9|9 8|5|5|7
## ind 9 8|1|6|2 8|2|5|5 8|6|6|3
## ind 10 1|6|6|4 2|5|4|4 1|9|2|3

4.3 Handling presence/absence data

adegenet was primarly designed to handle codominant, multiallelic markers like
microsatellites. However, dominant markers like AFLP can be used as well. In such a
case, only presence/absence of alleles can be deduced accurately from the genotypes. This
has several consequences, like the unability to compute allele frequencies. Hence, some
functionalities in adegenet won’t be available for dominant markers.
From version 1.2-3 of adegenet, the distinction between both types of markers is made by
the slot @type of genind or genpop objects, which equals codom for codominant markers, and
PA for presence/absence data. In the latter case, the ’tab’ slot of a genind object no longer
contains allele frequencies, but only presence/absence of alleles in a genotype. Similarly, the
tab slot of a genpop object not longer contains counts of alleles in the populations; instead,
it contains the number of genotypes in each population possessing at least one copy of the
concerned alleles. Moreover, in the case of presence/absence, the slots ’[Link]’, ’[Link]’,
and ’[Link]’ become useless, and are thus all set to NULL.

Objects of type ’PA’ are otherwise handled like usual (type ’codom’) objects. Operations
that are not available for PA type will issue an appropriate error message.
Here is an example using a toy dataset ’[Link]’ that can be downloaded from the
adegenet website, section ’Documentation’:

dat <- [Link]([Link]("files/[Link]",package="adegenet"), header=TRUE)

dat

## loc1 loc2 loc3 loc4

## indA 1 0 1 1

17
## indB 0 1 1 1
## indC 1 1 0 1
## indD 0 NA 1 NA
## indE 1 1 0 0
## indF 1 0 1 1
## indG 0 1 1 0

The function df2genind is used to obtain a genind object:

obj <- df2genind(dat, ploidy=1, type="PA")

obj

## /// GENIND OBJECT /////////

##
## // 7 individuals; 4 loci; 4 alleles; size: 4.2 Kb
##
## // Basic content
## @tab: 7 x 4 matrix of allele counts
## @[Link]: number of alleles per locus (range: 4-4)
## @ploidy: ploidy of each individual (range: 1-1)
## @type: PA
## @call: df2genind(X = dat, ploidy = 1, type = "PA")
##
## // Optional content
## - empty -

tab(obj)

## loc1 loc2 loc3 loc4

## indA 1 0 1 1
## indB 0 1 1 1
## indC 1 1 0 1
## indD 0 NA 1 NA
## indE 1 1 0 0
## indF 1 0 1 1
## indG 0 1 1 0

One can see that for instance, the summary of this object is more simple (no numbers of
alleles per locus, no heterozygosity):

pop(obj) <- rep(c('a','b'),4:3)

print(summary(obj))

## $n
## [1] 7

18
##
## $[Link]
## a b
## 4 3
##
## $[Link]
## [1] 7.142857

But we can still perform basic manipulation, like converting our object into a genpop:

obj2 <- genind2genpop(obj)

##
## Converting data from a genind to a genpop object...
##
## ...done.

obj2

## /// GENPOP OBJECT /////////

##
## // 2 populations; 4 loci; 4 alleles; size: 2.9 Kb
##
## // Basic content
## @tab: 2 x 4 matrix of allele counts
## @[Link]: number of alleles per locus (range: 4-4)
## @ploidy: ploidy of each individual (range: 1-1)
## @type: PA
## @call: genind2genpop(x = obj)
##
## // Optional content
## - empty -

tab(obj2)

## loc1 loc2 loc3 loc4

## a 2 2 3 3
## b 2 2 2 1

To continue with the toy example, we can perform a simple PCA. Allele presence absence
data are extracted and NAs replaced using tab:

X <- tab(obj, [Link]="mean")

Now the PCA is performed and plotted:

19
## make PCA
pca1 <- [Link](X,scannf=FALSE,scale=FALSE)
temp <- [Link](pop(obj))
myCol <- transp(c("blue","red"),.7)[temp]
myPch <- c(15,17)[temp]

## basic plot
plot(pca1$li, col=myCol, cex=3, pch=myPch)

## use wordcloud for non-overlapping labels

library(wordcloud)

## Error in library(wordcloud): there is no package called ’wordcloud’

textplot(pca1$li[,1], pca1$li[,2], words=rownames(X), cex=1.4, new=FALSE)

## Error in textplot(pca1$li[, 1], pca1$li[, 2], words = rownames(X), cex =

1.4, : could not find function "textplot"

## legend the axes by adding loadings

abline(h=0,v=0,col="grey",lty=2)
[Link](pca1$c1*.5, [Link]=TRUE)
legend("topright", pch=c(15,17), col=transp(c("blue","red"),.7),
leg=c("Group A","Group B"), [Link]=2)

20
 Group A

1.0
Group B
0.5
Axis2

loc1

loc2
0.0

loc4

loc3
−0.5

−0.5 0.0 0.5

Axis1

More generally, multivariate analyses from ade4, sPCA (spca), DAPC (dapc), the global
and local tests ([Link], [Link]), or the Monmonier’s algorithm (monmonier)
will work just fine with presence/absence data. However, it is clear that the usual Euclidean
distance (used in PCA and sPCA), as well as many other distances, is not as accurate
to measure genetic dissimilarity using presence/absence data as it is when using allele
frequencies. The reason for this is that in presence/absence data, a part of the information
is simply hidden. For instance, two individuals possessing the same allele will be considered
at the same distance, whether they possess one or more copies of the allele. This might be
especially problematic in organisms having a high degree of ploidy.

4.4 SNPs data

In adegenet, SNP data can be handled in two different ways. For relatively small datasets
(up to a few thousand SNPs) SNPs can be handled as usual codominant markers such as
microsatellites using genind objects. In the case of genome-wide SNP data (from hundreds
of thousands to millions of SNPs), genind objects are no longer efficient representation
of the data. In this case, we use genlight objects to store and handle information with
maximum efficiency and minimum memory requirements. See the tutorial genomics for
more information. Below, we introduce only the case of SNPs handled using genind objects.

21
 The most convenient way to convert SNPs into a genind is using df2genind, which is
described in the previous section. Let dat be an input matrix, as can be read into R using
[Link] or [Link], with genotypes in row and SNP loci in columns.

dat <- matrix(sample(c("a","t","g","c"), 15, replace=TRUE),nrow=3)

rownames(dat) <- paste("genot.", 1:3)
colnames(dat) <- 1:5
dat

## 1 2 3 4 5
## genot. 1 "c" "t" "c" "t" "c"
## genot. 2 "g" "a" "t" "c" "c"
## genot. 3 "t" "g" "t" "t" "a"

obj <- df2genind(dat, ploidy=1)

tab(obj)

## 1.c 1.g 1.t 2.t 2.a 2.g 3.c 3.t 4.t 4.c 5.c 5.a
## genot. 1 1 0 0 1 0 0 1 0 1 0 1 0
## genot. 2 0 1 0 0 1 0 0 1 0 1 1 0
## genot. 3 0 0 1 0 0 1 0 1 1 0 0 1

obj is a genind containing the SNPs information, which can be used for further analysis
in adegenet.

4.5 Extracting polymorphism from DNA sequences

This section only covers the cases of relatively small datasets which can be handled efficiently
using genind objects. For bigger (near full-genomes) datasets, SNPs can be extracted from
fasta files into a genlight object using fasta2genlight. See the tutorial on genomics for
more information.

DNA sequences can be read into R using ape’s [Link] (fasta and clustal formats), or
adegenet’s fasta2DNAbin for a RAM-greedy implementation (fasta alignments only). Other
options include reading data directly from GenBank using [Link], or from other
public databases using the seqinr package and transforming the alignment object into a
DNAbin using [Link]. Here, we illustrate this approach by re-using the example of
[Link]. A connection to the internet is required, as sequences are read directly
from a distant database.

library(ape)
ref <- c("U15717", "U15718", "U15719", "U15720",
"U15721", "U15722", "U15723", "U15724")
myDNA <- [Link](ref)
myDNA

22
## 8 DNA sequences in binary format stored in a list.
##
## All sequences of same length: 1045
##
## Labels:
## U15717
## U15718
## U15719
## U15720
## U15721
## U15722
## ...
##
## Base composition:
## a c g t
## 0.267 0.351 0.134 0.247
## (Total: 8.36 kb)

class(myDNA)

## [1] "DNAbin"

myDNA <- [Link](myDNA)

Polymorphism can be characterized using [Link] and [Link]: the first

plots SNP density along the alignment, the second tests whether these SNPs are randomly
distributed. For instance:

[Link](myDNA,codon=FALSE)

23
 Distribution of SNPs in the genome

0.0015
density

0.0010

0.0005

0.0000

0 250 500 750 1000

Nucleotide position

By default, the function differentiates nucleotide positions:

[Link](myDNA)

24
 Distribution of SNPs in the genome

0.006

0.004

Codon position
density

1
2
3

0.002

0.000

0 250 500 750 1000

Nucleotide position

In adegenet, only polymorphic loci are conserved to form a genind object. This
conversion is achieved by DNAbin2genind. This function allows one to specify a threshold for
polymorphism; for instance, one could retain only SNPs for which the second largest allele
frequency is greater than 1% (using the polyThres argument). This is achieved using:

obj <- DNAbin2genind(myDNA, polyThres=0.01)

obj

## /// GENIND OBJECT /////////

##
## // 8 individuals; 155 loci; 318 alleles; size: 93.5 Kb
##
## // Basic content
## @tab: 8 x 318 matrix of allele counts
## @[Link]: number of alleles per locus (range: 2-4)
## @[Link]: locus factor for the 318 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 1-1)
## @type: codom

25
## @call: DNAbin2genind(x = myDNA, polyThres = 0.01)
##
## // Optional content
## - empty -

Here, out of the 1,045 nucleotides of the sequences, 318 SNPs where extracted and stored
as a genind object. Positions of the SNPs are stored as names of the loci:

head(locNames(obj))

## [1] "11" "13" "26" "31" "34" "39"

4.6 Extracting polymorphism from proteic sequences

Alignments of proteic sequences can be exploited in adegenet in the same way as DNA
sequences (see section above). Alignments are scanned for polymorphic sites, and only
those are retained to form a genind object. Loci correspond to the position of the
residue in the alignment, and alleles correspond to the different amino-acids (AA). Aligned
proteic sequences are stored as objects of class alignment in the seqinr package [1]. See
?[Link] for a description of this class. The function extracting polymorphic sites
from alignment objects is alignment2genind.
Its use is fairly simple. It is here illustrated using a small dataset of aligned proteic
sequences:

library(seqinr)
[Link] <- [Link](file=[Link]("sequences/[Link]",
package="seqinr"), format = "mase")
[Link]

## $nb
## [1] 6
##
## $nam
## [1] "Langur" "Baboon" "Human" "Rat" "Cow" "Horse"
##
## $seq
## $seq[[1]]
## [1] "-kifercelartlkklgldgykgvslanwvclakwesgynteatnynpgdestdygifqinsrywcnngkpgavdachis
##
## $seq[[2]]
## [1] "-kifercelartlkrlgldgyrgislanwvclakwesdyntqatnynpgdqstdygifqinshywcndgkpgavnachis
##
## $seq[[3]]
## [1] "-kvfercelartlkrlgmdgyrgislanwmclakwesgyntratnynagdrstdygifqinsrywcndgkpgavnachls

26
##
## $seq[[4]]
## [1] "-ktyercefartlkrngmsgyygvsladwvclaqhesnyntqarnydpgdqstdygifqinsrywcndgkpraknacgip
##
## $seq[[5]]
## [1] "-kvfercelartlkklgldgykgvslanwlcltkwessyntkatnynpssestdygifqinskwwcndgkpnavdgchvs
##
## $seq[[6]]
## [1] "-kvfskcelahklkaqemdgfggyslanwvcmaeyesnfntrafngknangssdyglfqlnnkwwckdnkrsssnacnim
##
##
## $com
## [1] ";empty description\n" ";\n" ";\n" ";\n"
##
## attr(,"class")
## [1] "alignment"

x <- alignment2genind([Link])
x

## /// GENIND OBJECT /////////

##
## // 6 individuals; 82 loci; 212 alleles; size: 57.2 Kb
##
## // Basic content
## @tab: 6 x 212 matrix of allele counts
## @[Link]: number of alleles per locus (range: 2-5)
## @[Link]: locus factor for the 212 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 1-1)
## @type: codom
## @call: alignment2genind(x = [Link])
##
## // Optional content
## @other: a list containing: com

The six aligned protein sequences ([Link]) have been scanned for polymorphic sites,
and these have been extracted to form the genind object x. Note that several settings
such as the characters corresponding to missing values (i.e., gaps) and the polymorphism
threshold for a site to be retained can be specified through the function’s arguments (see
?alignment2genind).
The names of the loci directly provides the indices of polymorphic sites:

27
head(locNames(x))

## [1] "3" "4" "5" "6" "9" "11"

The table of polymorphic sites can be reconstructed easily by:

tabAA <- genind2df(x)

dim(tabAA)

## [1] 6 82

tabAA[, 1:20]

## 3 4 5 6 9 11 12 15 16 17 18 19 21 22 24 28 30 32 33 34
## Langur i f e r l r t k l g l d y k v n v l a k
## Baboon i f e r l r t r l g l d y r i n v l a k
## Human v f e r l r t r l g m d y r i n m l a k
## Rat t y e r f r t r n g m s y y v d v l a q
## Cow v f e r l r t k l g l d y k v n l l t k
## Horse v f s k l h k a q e m d f g y n v m a e

The global AA composition of the polymorphic sites is given by:

table(unlist(tabAA))

##
## a d e f g h i k l m n p q r s t v w y
## 35 38 16 9 33 13 27 28 31 8 44 10 26 47 36 20 42 6 23

Now that polymorphic sites have been converted into a genind object, simple distances
can be computed between the sequences. Note that adegenet does not implement specific
distances for protein sequences, we only use the simple Euclidean distance. Fancier protein
distances are implemented in R; see for instance [Link] in the seqinr package, and
[Link] in the phangorn package.

D <- dist(tab(x))
D

## Langur Baboon Human Rat Cow

## Baboon 5.291503
## Human 6.000000 5.291503
## Rat 8.717798 8.124038 8.602325
## Cow 7.874008 8.717798 8.944272 10.392305
## Horse 11.313708 11.313708 11.224972 11.224972 11.747340

28
 This matrix of distances is small enough for one to interprete the raw numbers. However,
it is also very straightforward to represent these distances as a tree or in a reduced space.
We first build a Neighbor-Joining tree using the ape package:

library(ape)
tre <- nj(D)
par(xpd=TRUE)
plot(tre, type="unrooted", edge.w=2)
edgelabels(tex=round(tre$[Link],1), bg=rgb(.8,.8,1,.8))

Baboon

Cow
Human
2.4
2.9
4.8
0.4
0.8
2.6 0.6
Langur

6.8
4.4

Horse

Rat

The best possible planar representation of these Euclidean distances is achieved by

Principal Coordinate Analyses (PCoA), which in this case will give identical results to PCA
of the original (centred, non-scaled) data:

pco1 <- [Link](D, scannf=FALSE,nf=2)

[Link](pco1$li*1.1, clab=0, pch="")
textplot(pco1$li[,1], pco1$li[,2], words=rownames(pco1$li),
cex=1.4, new=FALSE, xpd=TRUE)

29
## Error in textplot(pco1$li[, 1], pco1$li[, 2], words = rownames(pco1$li),
: could not find function "textplot"

title("Principal Coordinate Analysis\n-based on proteic distances-")

d=2
Principal Coordinate Analysis
−based on proteic distances−

4.7 Using genind/genpop constructors

genind or genpop objects are best obtained using the procedures described above. However,
one can also build new genind or genpop objects using the constructor new(). In this case,
the function must be given as main input an object corresponding to the @tab slot: a matrix
of integers with individuals in row and alleles in columns, with columns being named as
’[marker].[allele]’. Here is an example for genpop using a dataset from ade4 :

data(microsatt)
str(microsatt, [Link] = 1) # a plain list with a data frame and three vectors

## List of 4
## $ tab :'[Link]': 18 obs. of 112 variables:

30
## $ [Link] : chr [1:9] "INRA5" "INRA32" "INRA35" "INRA63" ...
## $ [Link] : int [1:9] 8 15 11 10 17 10 14 15 12
## $ [Link]: chr [1:112] "133" "137" "139" "141" ...

microsatt$tab[10:15,12:15]

## INRA32.168 INRA32.170 INRA32.174 INRA32.176

## Mtbeliard 0 0 0 1
## NDama 0 0 0 12
## Normand 1 0 0 2
## Parthenais 8 5 0 3
## Somba 0 0 0 20
## Vosgienne 2 0 0 0

microsatt$tab contains alleles counts per populations, and can therefore be used to make
a genpop object. Moreover, column names are set as required, and row names are unique. It
is therefore safe to convert these data into a genpop using the constructor:

toto <- new("genpop", tab=microsatt$tab)

toto

## /// GENPOP OBJECT /////////

##
## // 18 populations; 9 loci; 112 alleles; size: 30.4 Kb
##
## // Basic content
## @tab: 18 x 112 matrix of allele counts
## @[Link]: number of alleles per locus (range: 8-17)
## @[Link]: locus factor for the 112 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: .local(.Object = .Object, tab = ..1)
##
## // Optional content
## - empty -

summary(toto)

##
## // Number of populations: 18
## // Number of alleles per locus: 8 15 11 10 17 10 14 15 12
## // Number of alleles per group: 39 69 51 59 52 41 34 48 46 47 43 56 57 52 49 64 56 67
## // Percentage of missing data: 0 %

31
4.8 Exporting data
The genind class tends to become a standard in population genetics packages. As of
adegenet 2.0.0, export functions towards hierfstat have been removed, as the package now
uses genind objects as a native class. Similarly, export towards the package genetics have
been removed, as adegenet now relies on pegas for basic population genetics.

A generic way to export data is to produce a [Link] of genotypes coded by character

strings. This is done by genind2df:

obj <- genind2df(nancycats)

obj[1:5,1:5]

## pop fca8 fca23 fca43 fca45

## cat.1 P01 <NA> 136146 139139 116120
## cat.2 P01 <NA> 146146 139145 120126
## cat.3 P01 135143 136146 141141 116116
## cat.4 P01 133135 138138 139141 116126
## cat.5 P01 133135 140146 141145 126126

This function is flexible; for instance, one can separate alleles by any character string:

genind2df(nancycats,sep="|")[1:5,1:5]

## pop fca8 fca23 fca43 fca45

## cat.1 P01 <NA> 136|146 139|139 116|120
## cat.2 P01 <NA> 146|146 139|145 120|126
## cat.3 P01 135|143 136|146 141|141 116|116
## cat.4 P01 133|135 138|138 139|141 116|126
## cat.5 P01 133|135 140|146 141|145 126|126

Note that tabulations can be obtained as follows using ’\t’ character.

32
5 Basics of data analysis
5.1 Manipulating the data
Data manipulation is meant to be particularly flexible in adegenet. First, as genind and
genpop objects are basically formed by a data matrix (the @tab slot), it is natural to subset
these objects like it is done with a matrix. The [ operator does this, forming a new object
with the retained genotypes/populations and alleles:
data(microbov)
toto <- genind2genpop(microbov)

##
## Converting data from a genind to a genpop object...
##
## ...done.

toto

## /// GENPOP OBJECT /////////

##
## // 15 populations; 30 loci; 373 alleles; size: 96.3 Kb
##
## // Basic content
## @tab: 15 x 373 matrix of allele counts
## @[Link]: number of alleles per locus (range: 5-22)
## @[Link]: locus factor for the 373 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: genind2genpop(x = microbov)
##
## // Optional content
## @other: a list containing: coun breed spe

popNames(toto)

## [1] "Borgou" "Zebu" "Lagunaire" "NDama" "Somba"

## [10] "Charolais" "Gascon" "Limousin" "MaineAnjou" "Montbel

titi <- toto[1:3, ]

popNames(titi)

## [1] "Borgou" "Zebu" "Lagunaire"

The object toto has been subsetted, keeping only the first three populations. Of course, any
subsetting available for a matrix can be used with genind and genpop objects. In addition,

33
we can subset loci directly using indices or logicals, in which case they refer to the output of
locNames:

nAll(titi)

## INRA63 INRA5 ETH225 ILSTS5 HEL5 HEL1 INRA35 ETH152 INRA23 ETH10 HE
## 9 7 12 5 11 9 7 12 13 9
## ETH185 HAUT24 HAUT27 TGLA227 TGLA126 TGLA122 TGLA53 SPS115
## 16 13 12 15 8 22 21 9

tata <- titi[, loc=c(1, 3)]

tata

## /// GENPOP OBJECT /////////

##
## // 3 populations; 2 loci; 21 alleles; size: 18.9 Kb
##
## // Basic content
## @tab: 3 x 21 matrix of allele counts
## @[Link]: number of alleles per locus (range: 9-12)
## @[Link]: locus factor for the 21 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: .local(x = x, i = i, j = j, loc = ..1, drop = drop)
##
## // Optional content
## @other: a list containing: coun breed spe

nAll(tata)

## INRA63 ETH225
## 9 12

Alternatively, one can subset loci using their explicit name:

locNames(titi)

## [1] "INRA63" "INRA5" "ETH225" "ILSTS5" "HEL5" "HEL1" "INRA35" "ETH152"

## [18] "INRA37" "BM1818" "ILSTS6" "MM12" "CSRM60" "ETH185" "HAUT24" "HAUT27"

hel5 <- titi[, loc="HEL5"]

hel5

## /// GENPOP OBJECT /////////

##

34
## // 3 populations; 1 locus; 11 alleles; size: 17.2 Kb
##
## // Basic content
## @tab: 3 x 11 matrix of allele counts
## @[Link]: number of alleles per locus (range: 11-11)
## @[Link]: locus factor for the 11 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: .local(x = x, i = i, j = j, loc = "HEL5", drop = drop)
##
## // Optional content
## @other: a list containing: coun breed spe

locNames(hel5)

## [1] "HEL5"

When subsetting individuals/samples, some alleles may not be included in the subset
anymore. In case you want these alleles to be dropped, use the drop = TRUE argument.

data(nancycats)

nAll(nancycats[1:3, ])

## fca8 fca23 fca43 fca45 fca77 fca78 fca90 fca96 fca37

## 16 11 10 9 12 8 12 12 18

nAll(nancycats[1:3, , drop = TRUE])

## fca8 fca23 fca43 fca45 fca77 fca78 fca90 fca96 fca37

## 2 2 3 3 2 2 3 1 2

To simplify the task of separating data by marker systematically, the function seploc
can be used. It returns a list of objects (optionnaly, of data matrices), each corresponding
to a marker:

sepCats <- seploc(nancycats)

class(sepCats)

## [1] "list"

names(sepCats)

## [1] "fca8" "fca23" "fca43" "fca45" "fca77" "fca78" "fca90" "fca96" "fca37"

35
sepCats$fca45

## /// GENIND OBJECT /////////

##
## // 237 individuals; 1 locus; 9 alleles; size: 33.8 Kb
##
## // Basic content
## @tab: 237 x 9 matrix of allele counts
## @[Link]: number of alleles per locus (range: 9-9)
## @[Link]: locus factor for the 9 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: .local(x = x)
##
## // Optional content
## @pop: population of each individual (group size range: 9-23)
## @other: a list containing: xy

identical(tab(sepCats$fca45), tab(nancycats[,loc="fca45"]))

## [1] TRUE

The object sepCats$fca45 only contains data of the marker fca45.

Following the same idea, seppop allows one to separate genotypes in a genind object by
population. For instance, we can separate genotype of cattles in the dataset microbov by
breed:

data(microbov)
obj <- seppop(microbov)
class(obj)

## [1] "list"

names(obj)

## [1] "Borgou" "Zebu" "Lagunaire" "NDama" "Somba"

## [10] "Charolais" "Gascon" "Limousin" "MaineAnjou" "Montbel

obj$Borgou

## /// GENIND OBJECT /////////

##
## // 50 individuals; 30 loci; 373 alleles; size: 143.1 Kb
##

36
## // Basic content
## @tab: 50 x 373 matrix of allele counts
## @[Link]: number of alleles per locus (range: 5-22)
## @[Link]: locus factor for the 373 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: .local(x = x, i = i, j = j, pop = ..1, treatOther = ..2, quiet = ..3,
## drop = drop)
##
## // Optional content
## @pop: population of each individual (group size range: 50-50)
## @other: a list containing: coun breed spe

The returned object obj is a list of genind objects each containing genotypes of a given breed.

A last, rather vicious trick is to separate data by population and by marker. This is easy
using lapply; one can first separate population then markers, or the contrary. Here, we
separate markers inside each breed in obj:

obj <- lapply(obj,seploc)

names(obj)

## [1] "Borgou" "Zebu" "Lagunaire" "NDama" "Somba"

## [10] "Charolais" "Gascon" "Limousin" "MaineAnjou" "Montbel

class(obj$Borgou)

## [1] "list"

names(obj$Borgou)

## [1] "INRA63" "INRA5" "ETH225" "ILSTS5" "HEL5" "HEL1" "INRA35" "ETH152"

## [18] "INRA37" "BM1818" "ILSTS6" "MM12" "CSRM60" "ETH185" "HAUT24" "HAUT27"

obj$Borgou$INRA63

## /// GENIND OBJECT /////////

##
## // 50 individuals; 1 locus; 9 alleles; size: 13.8 Kb
##
## // Basic content
## @tab: 50 x 9 matrix of allele counts
## @[Link]: number of alleles per locus (range: 9-9)
## @[Link]: locus factor for the 9 columns of @tab

37
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: .local(x = x)
##
## // Optional content
## @pop: population of each individual (group size range: 50-50)
## @other: a list containing: coun breed spe

For instance, obj$Borgou$INRA63 contains genotypes of the breed Borgou for the marker
INRA63.

Lastly, one may want to pool genotypes in different datasets, but having the same markers,
into a single dataset. This is more than just merging the @tab components of all datasets,
because alleles can differ (they almost always do) and markers are not necessarily sorted
the same way. The function repool is designed to avoid these problems. It can merge any
genind provided as arguments as soon as the same markers are used. For instance, it can be
used after a seppop to retain only some populations:

obj <- seppop(microbov)

names(obj)
## [1] "Borgou" "Zebu" "Lagunaire" "NDama" "Somba"
## [10] "Charolais" "Gascon" "Limousin" "MaineAnjou" "Montbel
newObj <- repool(obj$Borgou, obj$Charolais)
newObj
## /// GENIND OBJECT /////////
##
## // 105 individuals; 30 loci; 295 alleles; size: 182.4 Kb
##
## // Basic content
## @tab: 105 x 295 matrix of allele counts
## @[Link]: number of alleles per locus (range: 4-17)
## @[Link]: locus factor for the 295 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: repool(obj$Borgou, obj$Charolais)
##
## // Optional content
## @pop: population of each individual (group size range: 50-55)
popNames(newObj)
## [1] "Borgou" "Charolais"

38
 Done !

Note that the content of @other can be processed during the conversion from genind to
genpop if the argument [Link] is set to TRUE. Only vectors of a length, or matrices
with a number of rows matching the number individuals will be processed. The way they
are processed is defined by a function passed as the argument [Link] (defaulting to
’mean’). Let us illustrate this using sim2pop:

data(sim2pop)
sim2pop

## /// GENIND OBJECT /////////

##
## // 130 individuals; 20 loci; 241 alleles; size: 192.1 Kb
##
## // Basic content
## @tab: 130 x 241 matrix of allele counts
## @[Link]: number of alleles per locus (range: 7-17)
## @[Link]: locus factor for the 241 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: old2new(object = sim2pop)
##
## // Optional content
## @pop: population of each individual (group size range: 30-100)
## @other: a list containing: xy

nInd(sim2pop)

## [1] 130

head(other(sim2pop)$xy)

## x y
## [1,] 35.11291 99.595997
## [2,] 22.57033 6.682107
## [3,] 76.99371 51.900514
## [4,] 44.31948 18.037868
## [5,] 94.40902 82.948821
## [6,] 51.29493 25.007193

dim(other(sim2pop)$xy)

## [1] 130 2

39
 The component sim2pop@other$xy contains spatial coordinates of individuals from 2
populations.

other(genind2genpop(sim2pop, [Link]=TRUE))

##
## Converting data from a genind to a genpop object...
##
## ...done.
## $xy
## x y
## P01 58.43405 48.37065
## P02 19.06501 61.00044

In this case, numeric vectors with a length corresponding to the number of individuals
will we averaged per groups; note that any other function than mean can be used by
providing any function to the argument [Link]. Matrices with a number of rows
corresponding to the number of individuals are processed similarly.

5.2 Using summaries

Both genind and genpop objects have a summary providing basic information about data.
Informations are both printed and invisibly returned as a list.

toto <- summary(nancycats)

names(toto)

## [1] "n" "[Link]" "[Link]" "[Link]" "[Link]" "Hobs" "Hexp"

par(mfrow=c(2,2))

plot(toto$[Link], toto$[Link], xlab="Colonies sample size",

ylab="Number of alleles",main="Alleles numbers and sample sizes",
type="n")
text(toto$[Link],toto$[Link],lab=names(toto$[Link]))

barplot(toto$[Link], ylab="Number of alleles",

main="Number of alleles per locus")

barplot(toto$Hexp-toto$Hobs, main="Heterozygosity: expected-observed",

ylab="Hexp - Hobs")

barplot(toto$[Link], main="Sample sizes per population",

ylab="Number of genotypes",las=3)

40
 Alleles numbers and sample sizes Number of alleles per locus

P11
P04

15
65
Number of alleles

Number of alleles
P14
P06
55

10
P08 P02
P12
P03
P05
P10
45

5
P09 P15 P13
P07
P16
P01
35

P17

0
10 12 14 16 18 20 22 fca8 fca43 fca77 fca90 fca37

Colonies sample size

Heterozygosity: expected−observed Sample sizes per population

0.20

20
Number of genotypes
0.15
Hexp − Hobs

15
0.10

10
0.05

5
0.00

P01
P02
P03
P04
P05
P06
P07
P08
P09
P10
P11
P12
P13
P14
P15
P16
P17
fca8 fca43 fca77 fca90 fca37

par(mfrow = c(1,1))

Is mean observed H significantly lower than mean expected H ?

[Link](list(toto$Hexp,toto$Hobs))

##
## Bartlett test of homogeneity of variances
##
## data: list(toto$Hexp, toto$Hobs)
## Bartlett's K-squared = 0.046962, df = 1, p-value = 0.8284

[Link](toto$Hexp,toto$Hobs,pair=T,[Link]=TRUE,alter="greater")

##
## Paired t-test
##
## data: toto$Hexp and toto$Hobs

41
## t = 8.3294, df = 8, p-value = 1.631e-05
## alternative hypothesis: true mean difference is greater than 0
## 95 percent confidence interval:
## 0.1134779 Inf
## sample estimates:
## mean difference
## 0.1460936

Yes, it is.

5.3 Testing for Hardy-Weinberg equilibrium

As of version 2.0.0, adegenet is designed to work alongside a number of other packages,
especially pegas and hierfstat for a number of classical population genetics methods. The
former function [Link] has consequently been removed, and replaced by pegas’s
[Link], which performs one test per locus:

library(pegas)
data(nancycats)
[Link] <- [Link](nancycats, B=0)
[Link]

## chi^2 df Pr(chi^2 >)

## fca8 395.80006 120 0.000000e+00
## fca23 239.34221 55 0.000000e+00
## fca43 434.33397 45 0.000000e+00
## fca45 66.11849 36 1.622163e-03
## fca77 270.52066 66 0.000000e+00
## fca78 402.80002 28 0.000000e+00
## fca90 217.19836 66 0.000000e+00
## fca96 193.36764 66 1.965095e-14
## fca37 291.00731 153 1.209777e-10

Note that B=0 is used for the parametric version; larger numbers will indicate the number
of permutations to use for a Monte-Carlo version.

5.4 Measuring and testing population structure (a.k.a F statistics)

Population structure is traditionally measured and tested using F statistics, in particular
Fst . As of version 2.0.0, adegenet relies on hierfstat and pegas for most F statistics.
Can we find any population structure in the cat colonies from Nancy? The basic Weir
and Cockerham [12] F statistics are provided by the wc() function from heirfstat:

42
library("hierfstat")

##
## Attaching package: ’hierfstat’
## The following objects are masked from ’package:adegenet’:
##
## Hs, [Link]
## The following objects are masked from ’package:ape’:
##
## pcoa, varcomp

wc(nancycats)

## $FST
## [1] 0.08494959
##
## $FIS
## [1] 0.120589

This table provides two F statistics F st (pop/total), and F is (ind/pop). These are
overall measures which take into account all genotypes and all loci.

For more detail, pegas provides estimates by locus:

library(pegas)
ftab <- Fst([Link](nancycats))
ftab # per-locus F-statistics

## Fit Fst Fis

## fca8 0.23508749 0.10150515 0.148673460
## fca23 0.16462947 0.06746762 0.104191391
## fca43 0.15144873 0.06893755 0.088620458
## fca45 0.07518537 0.07652596 -0.001451681
## fca77 0.27904947 0.10036588 0.198618075
## fca78 0.18424901 0.07025915 0.122603911
## fca90 0.20987443 0.09168833 0.130116240
## fca96 0.19425217 0.10981110 0.094857474
## fca37 0.26040330 0.06985321 0.204860244

colMeans(ftab) # global F-statistics

## Fit Fst Fis

## 0.19490883 0.08404599 0.12123217

Confidence intervals for these F statistics can be obtained through the [Link]() function
in hierfstat, which takes a data frame of population strata and a data frame of genotypes.

43
You can convert the genind object to this data frame with genind2hierfstat()

nc <- genind2hierfstat(nancycats)
[Link](nc[1], nc[-1])$ci

## H-Total F-pop/Total F-Ind/Total H-pop F-Ind/pop Hobs

## 2.5% 0.7276 0.0751 0.1541 0.6683 0.0802 0.5754
## 50% 0.7814 0.0852 0.1970 0.7147 0.1210 0.6296
## 97.5% 0.8294 0.0962 0.2339 0.7538 0.1576 0.6679

Finally, pairwise Fst is frequently used as a measure of distance between populations. The
function [Link] can compute Nei’s estimator [10] of pairwise F st, defined as:

Ht − (nA Hs (A) + nB Hs (B))/(nA + nB )

F st(A, B) =
Ht
where A and B refer to the two populations of sample size nA and nB and respective expected
heterozygosity Hs (A) and Hs (B), and Ht is the expected heterozygosity
P 2 in the whole dataset.
For a given locus, expectedPheterozygosity is computed as 1 − p i , where pi is the frequency
of the ith allele, and the represents summation over all alleles. For multilocus data, the
heterozygosity is simply averaged over all loci. These computations are achieved for all pairs
of populations by the function [Link]() with the “Nei87” method; we illustrate this on
a subset of individuals of nancycats (computations for the whole dataset would take a few
tens of seconds):

matFst <- [Link](nancycats[1:50, ], method = "Nei87")

matFst

## P01 P02 P03

## P02 0.1334
## P03 0.0855 0.1307
## P04 0.0831 0.1174 0.0166

The resulting matrix is Euclidean when there are no missing values:

[Link](matFst)

## [1] TRUE

It can therefore be used in a Principal Coordinate Analysis (which requires Euclideanity),

used to build trees, etc.

5.5 Estimating inbreeding

Inbreeding refers to an excess of homozygosity in a given individual due to the mating of
genetically related parents. This excess of homozygosity is due to the fact that there are

44
non-negligible chances of inheriting two identical alleles from a recent common ancestor.
Inbreeding can be associated to a loss of fitness leading to ”inbreeding depression”. Typically,
loss of fitness is caused by recessive deleterious alleles which have usually low frequency in
the population, but for which inbred individuals are more likely to be homozygotes.

The inbreeding coefficient F is defined as the probability that at a given locus, two
identical alleles have been inherited from a common ancestor. In the absence ofP inbreeding,
the probability of being homozygote at one loci is (for diploid individuals) simply i p2i where
i indexes the alleles and pi is the frequency of allele i. This can be generalized incorporating
F as: X
p(homozygote) = F + (1 − F ) p2i
i

and even more generally, for any ploidy π:

X
p(homozygote) = F + (1 − F ) pπi
i

This therefore allows for computing the likelihood of a given state (homozygote/heterozygote)
in a given genotype (log-likelihood are summed across loci for more than one marker).

This estimation is achieved by inbreeding. Depending on the value of the argument

[Link], the function returns a sample from the likelihood function ([Link]=’sample’)
or the likelihood function itself, as a R function ([Link]=’function’). While likelihood
functions are quickly obtained and easy to display graphically, sampling from the distributions
is more computer intensive but useful to derive summary statistics of the distributions. Here,
we illustrate inbreeding using the microbov dataset, which contains cattle breeds genotypes
for 30 microsatellites; to focus on breed Salers only, we use seppop:

data(microbov)
sal <- seppop(microbov)$Salers
sal

## /// GENIND OBJECT /////////

##
## // 50 individuals; 30 loci; 373 alleles; size: 143.1 Kb
##
## // Basic content
## @tab: 50 x 373 matrix of allele counts
## @[Link]: number of alleles per locus (range: 5-22)
## @[Link]: locus factor for the 373 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: .local(x = x, i = i, j = j, pop = ..1, treatOther = ..2, quiet = ..3,
## drop = drop)

45
##
## // Optional content
## @pop: population of each individual (group size range: 50-50)
## @other: a list containing: coun breed spe

We first compute the mean inbreeding for each individual, and plot the resulting
distribution:

temp <- inbreeding(sal, N=100)

class(temp)

## [1] "list"

head(names(temp))

## [1] "FRBTSAL9087" "FRBTSAL9088" "FRBTSAL9089" "FRBTSAL9090" "FRBTSAL9091" "FRBTSAL909

head(temp[[1]],20)

## [1] 0.08652402 0.11774134 0.16216700 0.11806648 0.15597850 0.09334167 0.14755845 0.2

## [17] 0.05989588 0.10744951 0.05727140 0.03256880

temp is a list of values sampled from the likelihood distribution of each individual; means
values are obtained for all individuals using sapply:

Fbar <- sapply(temp, mean)

hist(Fbar, col="firebrick", main="Average inbreeding in Salers cattles")

46
 Average inbreeding in Salers cattles

15
10
Frequency

5
0

0.1 0.2 0.3 0.4

Fbar

We can see that some individuals (actually, a single one) have higher inbreeding (>0.4).
We can recompute inbreeding for this individual, asking for the likelihood function to be
returned:

which(Fbar>0.4)

## FRBTSAL9266
## 37

F <- inbreeding(sal, [Link]="function")[which(Fbar>0.4)]

F

## $FRBTSAL9266
## function (Fest)
## {
## args <- lapply([Link]([Link]())[-1L], eval, [Link]())
## names <- if ([Link](names(args)))
## character(length(args))
## else names(args)

47
## dovec <- names %in% [Link]
## [Link]("mapply", c(FUN = FUN, args[dovec], MoreArgs = list(args[!dovec]),
## SIMPLIFY = SIMPLIFY, [Link] = [Link]))
## }
## <environment: 0x7fd95adeb628>

The output object F can seem a bit cryptic: it is an function embedded within a hidden
environment. This does not matter, however, since it is easily represented:

plot(F$FRBTSAL9266, main=paste("Inbreeding of individual",names(F)),

xlab="Inbreeding (F)", ylab="Probability density")

Inbreeding of individual FRBTSAL9266

8e−09
6e−09
Probability density

4e−09
2e−09
0e+00

0.0 0.2 0.4 0.6 0.8 1.0

Inbreeding (F)

Indeed, this individual shows subsequent inbreeding, with about 50% chances of being
homozygote through inheritance from a common ancestor of its parents.

48
6 Multivariate analysis
6.1 General overview
Multivariate analysis consists in summarising a strongly multivariate information into a few
synthetic variables. In genetics, such approaches are useful to get a simplified picture of
the genetic diversity obersved amongst individuals or populations. A review of multivariate
analysis in population genetics can be found in [6]. Here, we aim at providing an overview
of some applications using methods implemented in ade4 and adegenet.

Useful functions include:

ˆ tab (adegenet): extract allele counts or frequencies and replaces missing data; useful,
among other things, before running a principal component analysis (PCA).

ˆ [Link] (ade4 ): implements PCA; can be used on transformed allele frequencies of

individuals or populations.

ˆ [Link] (ade4 ): implements Correspondance Analysis (CA); can be used on raw allele
counts of populations (@tab slot in genpop objects).

ˆ [Link] (adegenet): implements 5 pairwise genetic distances between populations

ˆ [Link] (hierfstat): implements pairwise FST , which is also a Euclidean distance

between populations.

ˆ dist (stats): computes pairwise distances between multivariate observations; can be

used on raw or transformed allele frequencies.

ˆ [Link] (ade4 ): implements Principal Coordinates Analysis (PCoA); this methods

finds synthetic variables which summarize a Euclidean distance matrix as best as
possible; can be used on outputs of dist, [Link], and [Link].

ˆ [Link] (ade4 ): tests whether a distance matrix is Euclidean, which is a pre-requisite

of PCoA.

ˆ cailliez (ade4 ): renders a non-Euclidean distance matrix Euclidean by adding a

constant to all entries.

ˆ dapc (adegenet): implements the Discriminant Analysis of Principal Components

(DAPC [7]), a powerful method for the analysis of population genetic structures; see
dedicated tutorial (dapc).

ˆ sPCA (adegenet): implements the spatial Principal Component Analysis (sPCA [4]), a
method for the analysis of spatial genetic structures; see dedicated tutorial (dapc).

ˆ glPca (adegenet): implements PCA for genome-wide SNP data stored as genlight
objects; see dedicated tutorial (genomics).

49
 Besides the procedures themselves, graphic functions are also often of the utmost
importance; these include:

ˆ scatter (ade4,adegenet): generic function to display multivariate analyses; in practice,

the most useful application for genetic data is the one implemented in adegenet for
DAPC results.

ˆ [Link] (ade4 ): function used for basic display of principal components.

ˆ loadingplot (adegenet): function used to display the loadings (i.e., contribution to a

given structure) of alleles for a given principal component; annotates and returns the
most contributing alleles.

ˆ [Link] (ade4 ): displays two quantitative variables with known groups of observations,
using inertia ellipses for the groups; useful to represent principal components when
groups are known.

ˆ [Link] (ade4 ): same as [Link], except convex polygons are used rather than ellipses.

ˆ [Link] (ade4 ): graphical display of a quantitative variable distributed over a two-

dimensional space; useful to map principal components or allele frequencies over a
geographic area.

ˆ colorplot (adegenet): graphical display of 1 to 3 quantitative variables distributed over

a two-dimensional space; useful for combined representations of principal components
over a geographic area. Can also be used to produce color versions of traditional
scatterplots.

ˆ transp (adegenet): auxiliary function making colors transparent.

ˆ num2col (adegenet): auxiliary function transforming a quantitative variable into colors

using a given palette.

ˆ assignplot (adegenet): specific plot of group membership probabilities for DAPC; see
dedicated tutorial (dapc).

ˆ compoplot (adegenet): specific ’STRUCTURE-like’ plot of group membership

probabilities for DAPC; see dedicated tutorial (dapc).

ˆ [Link] (ade4 ): add inset plots to an existing figure.

ˆ [Link] (ade4 ): specific application of [Link] to add barplots of

eigenvalues to an existing figure.

In the sections below, we briefly illustrate how these tools can be combined to extract
information from genetic data.

50
6.2 Performing a Principal Component Analysis on genind objects
The tables contained in genind objects can be submitted to a Principal Component Analysis
(PCA) to seek a summary of the genetic diversity among the sampled individuals. Such
analysis is straightforward using adegenet to prepare data and ade4 for the analysis per
se. One has first to extract allele counts or frequencies from the genind object and replace
missing data (NAs) by the mean allele frequency. This is achieved by tab:

data(microbov)
sum([Link](microbov$tab))

## [1] 6325

There are 6325 missing data, which will be replaced by tab:

X <- tab(microbov, freq = TRUE, [Link] = "mean")

class(X)

## [1] "matrix" "array"

dim(X)

## [1] 704 373

X[1:5,1:5]

## INRA63.167 INRA63.171 INRA63.173 INRA63.175 INRA63.177

## AFBIBOR9503 0 0 0 0 0.0
## AFBIBOR9504 0 0 0 0 0.0
## AFBIBOR9505 0 0 0 0 0.5
## AFBIBOR9506 0 0 0 0 0.0
## AFBIBOR9507 0 0 0 0 0.5

The analysis can now be performed. We disable the scaling in [Link], as all ’variables’
(alleles) are vary on a common scale. Note: in practice, retained axes can be chosen
interactively by removing the arguments scannf=FALSE,nf=3.

pca1 <- [Link](X, scale = FALSE, scannf = FALSE, nf = 3)

barplot(pca1$eig[1:50], main = "PCA eigenvalues", col = [Link](50))

51
 PCA eigenvalues

1.2
1.0
0.8
0.6
0.4
0.2
0.0

pca1

## Duality diagramm
## class: pca dudi
## $call: [Link](df = X, scale = FALSE, scannf = FALSE, nf = 3)
##
## $nf: 3 axis-components saved
## $rank: 341
## eigen values: 1.27 0.5317 0.423 0.2853 0.2565 ...
## vector length mode content
## 1 $cw 373 numeric column weights
## 2 $lw 704 numeric row weights
## 3 $eig 341 numeric eigen values
##
## [Link] nrow ncol content
## 1 $tab 704 373 modified array
## 2 $li 704 3 row coordinates
## 3 $l1 704 3 row normed scores

52
## 4 $co 373 3 column coordinates
## 5 $c1 373 3 column normed scores
## other elements: cent norm

The output object pca1 is a list containing various information; of particular interest are:
ˆ $eig: the eigenvalues of the analysis, indicating the amount of variance represented by
each principal component (PC).
ˆ $li: the principal components of the analysis; these are the synthetic variables
summarizing the genetic diversity, usually visualized using scatterplots.
ˆ $c1: the allele loadings, used to compute linear combinations forming the PCs; squared,
they represent the contribution to each PCs.
The basic scatterplot for this analysis can be obtained by:

[Link](pca1$li)
title("PCA of microbov dataset\naxes 1-2")
[Link](pca1$eig[1:20], 3,1,2)
d=1
PCA of microbov dataset
axes 1−2 AFBIZEB9499
AFBIZEB9477
AFBIZEB9454
AFBIZEB9494 AFBIZEB9453
AFBIZEB9502
AFBIZEB9474
AFBIZEB9455
AFBIZEB9488
AFBIZEB9459AFBIZEB9493
AFBIZEB9501
AFBIZEB9496
AFBIZEB9456
AFBIZEB9484
AFBIZEB9469 AFBIBOR9546
AFBIZEB9498
AFBIZEB9479
AFBIBOR9550
AFBIZEB9473
AFBIZEB9492
AFBIZEB9486
AFBIZEB9460
AFBIZEB9500
AFBIZEB9490
AFBIZEB9457
AFBIZEB9462
AFBIZEB9475
AFBIBOR9503
AFBIZEB9487
AFBIZEB9480
AFBIZEB9470
AFBIZEB9464
AFBIZEB9482 AFBIZEB9497
AFBIZEB9472
AFBIZEB9478
AFBIZEB9483
AFBIZEB9485
AFBIZEB9466
AFBIZEB9471 AFBIZEB9495
AFBIBOR9544
AFBIZEB9467 AFBIZEB9481
AFBIZEB9491
AFBIZEB9468
AFBIBOR9537
AFBIZEB9476
AFBIBOR9520
AFBIBOR9551
AFBIBOR9504AFBIZEB9489
AFBIBOR9534
AFBIBOR9547
AFBIZEB9461
AFBIBOR9552
AFBIBOR9525
AFBIBOR9524
AFBIBOR9526
AFBIZEB9463 AFBIBOR9545
FRBTAUB9071 AFBIZEB9458
AFBIBOR9522
AFBIBOR9518 AFBIBOR9540
AFBIBOR9538
AFBIBOR9536
AFBIBOR9535
FRBTGAS9173 AFBIBOR9515
AFBIBOR9514 AFBIBOR9539
AFBIBOR9506
AFBIBOR9532AFBIBOR9513
FRBTSAL9258 AFBIBOR9519 AFBIBOR9530
AFBIZEB9465 AFBIBOR9527
AFBIBOR9543AFBIBOR9529
AFBIBOR9531
AFBIBOR9533
FRBTAUB9289
FRBTSAL9102 FRBTBDA29873 FRBTMA25382
FRBTMBE1544
FRBTAUB9072 AFBIBOR9516AFBIBOR9523
AFBIBOR9549
AFBIBOR9507
AFBIBOR9512
FRBTAUB9224 FRBTSAL9245
FRBTAUB9063
FRBTSAL9100
FRBTAUB9218
FRBTSAL9241 FRBTBPN1875FRBTLIM5136 AFBIBOR9528
AFBIBOR9521
AFBIBOR9542
AFBIBOR9517
FRBTAUB9223
FRBTGAS9054 FRBTAUB9062
FRBTLIM30855
FRBTAUB9208
FRBTBDA29874
FRBTAUB9238
FRBTGAS14186FRBTSAL9259
FRBTAUB9237
FRBTAUB9225
FRBTSAL9251
FRBTSAL9257 FRBTGAS9203
FRBTAUB9066
FRBTSAL9267
FRBTMA29809
FRBTSAL9087 FRBTGAS9184
FRBTAUB9061 FRBTGAS9199
FRBTBPN1902
FRBTCHA26797 AFBIBOR9548
FRBTGAS9200
FRBTSAL9265
FRBTLIM30852
FRBTSAL9088
FRBTAUB9234
FRBTLIM30837
FRBTGAS9201 FRBTBDA29875
FRBTCHA30341
FRBTBDA35931
FRBTCHA26193
FRBTBDA35280
FRBTAUB9211
FRBTBDA35244 FRBTBDA29870
FRBTGAS9202
FRBTCHA26196 AFBTSOM9378 AFBTSOM9386
FRBTSAL9243 FRBTLIM30856
FRBTMA25418
FRBTSAL9091
FRBTSAL9099
FRBTAUB9221
FRBTBPN1935
FRBTAUB9212 FRBTBDA35485
FRBTBDA35258
FRBTAUB9081
FRBTSAL9271
FRBTBPN1901
FRBTBDA29878 FRBTGAS9188
FRBTSAL9097
FRBTAUB9065 FRBTBDA35877
FRBTGAS9171
FRBTLIM30823
FRBTCHA30879
FRBTBPN25811
FRBTGAS9195
FRBTBDA35446
FRBTCHA26202
FRBTAUB9290
FRBTGAS9180
FRBTLIM30838
FRBTBPN1877 AFBTND208AFBTSOM9393
AFBIBOR9541
AFBTSOM9371
AFBTSOM9383
FRBTSAL9262
FRBTLIM30839FRBTSAL9255
FRBTSAL9095
FRBTSAL9270 FRBTGAS9049
FRBTBDA29861 FRBTCHA30391
FRBTGAS9056
FRBTBDA29855
FRBTMBE1523
FRBTBDA35864
FRBTBDA35379
FRBTCHA26783
FRBTCHA25015
FRBTBDA29865
FRBTLIM30857
FRBTGAS14188
FRBTCHA26054
FRBTBPN1927
FRBTMA25278
FRBTAUB9077
FRBTAUB9069
FRBTBDA35284
FRBTAUB9229
FRBTBPN1906
FRBTGAS9183 FRBTCHA30335 AFBIBOR9509
FRBTAUB9219
FRBTSAL9283
FRBTCHA25513
FRBTLIM30854
FRBTBDA36124 FRBTGAS9179
FRBTAUB9073
FRBTGAS14183
FRBTGAS9060
FRBTBDA35286
FRBTSAL9277
FRBTSAL9252
FRBTGAS9055
FRBTAUB9232
FRBTSAL9285
FRBTBDA29881 FRBTMA25820
FRBTMA29948
FRBTLIM5135
FRBTCHA30382
FRBTBDA29853
FRBTLIM30834
FRBTBDA35274
FRBTGAS9174
FRBTBDA35256
FRBTCHA26785
FRBTBDA35916
FRBTBDA35259 FRBTCHA25295
FRBTBDA29857
FRBTCHA26274
FRBTLIM30836
FRBTBDA29854 AFBIBOR9511
AFBIBOR9505
AFBTND241
FRBTSAL9275
FRBTSAL9096 FRBTGAS9186
FRBTLIM30851
FRBTBDA29864
FRBTAUB9074
FRBTAUB9216
FRBTSAL9261 FRBTMA25298
FRBTMA25594
FRBTMBE1517
FRBTLIM30827
FRBTAUB9230 FRBTMBE1505
FRBTLIM30840
FRBTAUB9231
FRBTAUB9227
FRBTMA29571
FRBTSAL9094
FRBTGAS14180
FRBTMA25412
FRBTMA25428
FRBTCHA26014
FRBTLIM30860
FRBTBAZ15655
FRBTAUB9075 FRBTGAS9050
FRBTAUB9068
FRBTAUB9070 FRBTMA25273
FRBTBPN1872
FRBTLIM30830
FRBTSAL9093
FRBTBPN1873
FRBTAUB9228 AFBIBOR9510 AFBTSOM9375
FRBTGAS9198
FRBTAUB9286
FRBTSAL9103 FRBTBAZ30436
FRBTLIM30829
FRBTSAL9276
FRBTSAL9246
FRBTBDA35262
FRBTSAL9260
FRBTGAS9176
FRBTSAL9101
FRBTSAL9256 FRBTMBE1536
FRBTCHA26205
FRBTGAS9181
FRBTCHA26786
FRBTAUB9220 FRBTMA25902
FRBTAUB9222 FRBTLIM30847
FRBTCHA15985
FRBTLIM30826
FRBTGAS9182 FRBTCHA25390
FRBTBDA35283
FRBTCHA26790 FRBTBDA29882
FRBTGAS9205
FRBTMBE1519
FRBTBDA29856
FRBTMA26019
FRBTCHA26793
FRBTBDA35245
FRBTBPN1911
FRBTBPN1897 FRBTCHA26789
FRBTBPN1894
FRBTMBE1506
FRBTCHA26792
FRBTLIM30819 FRBTBDA35243
FRBTSAL9284 FRBTMA26022
FRBTCHA26800
FRBTBPN1928 AFBTND212
AFBTSOM9359
FRBTAUB9226 FRBTBAZ26396
FRBTCHA26040
FRBTLIM30822
FRBTGAS9051
FRBTGAS9053
FRBTBDA29851
FRBTAUB9064
FRBTBAZ30539
FRBTAUB9067
FRBTSAL9273 FRBTMA29561
FRBTCHA30884
FRBTBPN1899
FRBTAUB9078
FRBTLIM30853
FRBTGAS9185
FRBTMA29819
FRBTGAS9172
FRBTLIM30820
FRBTLIM30843FRBTCHA26074
FRBTCHA25995
FRBTCHA30896
FRBTMA25612
FRBTCHA30373
FRBTBDA35269
FRBTLIM30828
FRBTMBE1511
FRBTBAZ25956 FRBTCHA25009
FRBTCHA25069
FRBTCHA26784
FRBTBDA29879 FRBTMA29812
FRBTGAS9189 FRBTBDA36120
FRBTBPN1910
FRBTBDA35899
FRBTLIM30859 AFBTSOM9357 AFBTSOM9401
AFBTND285
AFBIBOR9508
AFBTND217
AFBTSOM9368
FRBTSAL9247
FRBTMBE1548
FRBTLIM30835
FRBTBAZ26403 FRBTMBE1534
FRBTBAZ26463
FRBTLIM30858
FRBTGAS9052
FRBTLIM3001 FRBTBDA35268
FRBTMA25488
FRBTAUB9079
FRBTBPN1934 FRBTGAS9193
FRBTMA26035
FRBTCHA26285
FRBTBAZ29254
FRBTMA25982
FRBTMA26491
FRBTGAS9059
FRBTBDA35941
FRBTMA25439
FRBTMA29943
FRBTLIM30831
FRBTCHA15957
FRBTCHA26011
FRBTCHA25322
FRBTMBE1541 FRBTMBE1520
FRBTGAS14187
FRBTBDA29869 FRBTBDA29872
FRBTBDA29859
FRBTMBE1502
FRBTBAZ25576
FRBTSAL9280
FRBTLIM30825
FRBTBDA29858 FRBTMA26321
FRBTBDA35255
FRBTBDA29868
FRBTLIM30824
FRBTBAZ26375 FRBTGAS9057 AFBTSOM9369 AFBTSOM9390
FRBTSAL9242
FRBTMA25423
FRBTAUB9076
FRBTBAZ30533
FRBTAUB9235
FRBTAUB9288 FRBTSAL9272
FRBTBDA29880
FRBTLIM30848
FRBTBAZ26097
FRBTSAL9090 FRBTMA26232
FRBTMBE1508
FRBTCHA30886
FRBTBDA35267
FRBTBAZ26400 FRBTGAS9177
FRBTBDA35248
FRBTLIM30845
FRBTBAZ29262
FRBTMA25387
FRBTMA25978
FRBTCHA30344 FRBTLIM30818
FRBTCHA25654
FRBTLIM30844
FRBTBAZ29247
FRBTMBE1513
FRBTMBE1547
FRBTLIM30849
FRBTBDA35955
FRBTSAL9249
FRBTSAL9253
FRBTBDA35260 FRBTMBE1497
FRBTMA25436
FRBTBPN1870
FRBTLIM30816
FRBTBAZ26386
FRBTBAZ29244
FRBTMA25896
FRBTCHA30353 AFBTND233 AFBTSOM9354
AFBTSOM9396
FRBTCHA25543 FRBTBPN1904
FRBTMBE1540
FRBTBAZ30421
FRBTAUB9287
FRBTGAS9204
FRBTBAZ26457
FRBTBPN1895
FRBTBAZ30531
FRBTBPN1930 FRBTBPN1896
FRBTMA25922
FRBTMBE1510
FRBTCHA25326
FRBTBPN25810
FRBTLIM30842
FRBTMA25917
FRBTLIM5137
FRBTMA25530 AFBTND259
AFBTND292 AFBTSOM9394
AFBTSOM9391
AFBTSOM9374
AFBTSOM9360
AFBTSOM9365
AFBTND284
FRBTLIM30817
FRBTSAL9248
FRBTBAZ29272
FRBTBDA29862 FRBTBAZ15654
FRBTBAZ25578 FRBTMBE1531
FRBTBAZ25950
FRBTAUB9210
FRBTMBE1529
FRBTCHA25024
FRBTBAZ26112
FRBTGAS9190
FRBTMBE1496
FRBTGAS9197 FRBTLIM30846
FRBTBAZ26110
FRBTMBE1549
FRBTGAS14185
FRBTCHA25018
FRBTCHA30356 AFBTND209
AFBTSOM9364
AFBTSOM9363AFBTSOM9370
AFBTSOM9382
AFBTND242
FRBTBAZ25954 FRBTCHA15946
FRBTBPN1937
FRBTBPN1876
FRBTMA26280
FRBTMBE1535
FRBTGAS9178
FRBTBAZ29261
FRBTCHA26199
FRBTLIM30841
FRBTMA25684
FRBTLIM30832 FRBTBAZ26388 FRBTBDA29852
FRBTGAS9170
FRBTCHA30878
FRBTBPN1913
FRBTBAZ29598
FRBTSAL9089
FRBTBPN1914 AFBTND205
AFBTSOM9392AFBTSOM9356
FRBTMA29950
FRBTMA26168
FRBTMA29806
FRBTSAL9268
FRBTBAZ26092
FRBTMBE1538 FRBTBAZ30429
FRBTBDA29863
FRBTBPN1915
FRBTCHA30893
FRBTBAZ29281 FRBTBDA35281
FRBTBAZ30540
FRBTMA25522
FRBTBAZ26352
FRBTSAL9266 FRBTBAZ29253
FRBTCHA26798
FRBTLIM30833
FRBTGAS9058 AFBTSOM9380
AFBTSOM9395 AFBTSOM9387
AFBTSOM9366
AFBTSOM9373
AFBTND206
AFBTND258
AFBTND213
FRBTLIM5133 FRBTCHA26246
FRBTSAL9250
FRBTBAZ26354 FRBTBDA29860
FRBTLIM30850
FRBTMA25433
FRBTMA29572
FRBTMA29945
FRBTMA29952
FRBTBDA35278
FRBTGAS14184 FRBTBAZ26078
FRBTBDA29877
FRBTMA25588 AFBTSOM9355
AFBTSOM9398
AFBTND202
AFBTSOM9381
FRBTMA25614
FRBTCHA15994
FRBTBPN1932
FRBTBAZ30440
FRBTGAS9175
FRBTMBE1516 AFBTND214
AFBTND223 AFBTSOM9389
AFBTND254
FRBTBAZ26439 FRBTMA25409
FRBTMBE1507
FRBTBAZ25957 FRBTBPN1898
FRBTBPN1907
FRBTBDA29866 FRBTAUB9213
FRBTBAZ29259
FRBTLIM30821 AFBTSOM9367
AFBTSOM9379 AFBTND207
AFBTSOM9397
AFBTSOM9353
AFBTSOM9384
FRBTBAZ26469 AFBTSOM9399
AFBTSOM9361
AFBTSOM9376 AFBTSOM9377
AFBTSOM9400
AFBTSOM9388
FRBTMBE1532FRBTMBE1514
FRBTBAZ30420
FRBTBAZ26401
FRBTBAZ30418 AFBTSOM9362AFBTND248
AFBTND216
AFBTSOM9372
AFBTND257 AFBTND253
AFBTLAG9409
AFBTSOM9352
FRBTMA25282 FRBTBAZ29246
FRBTMBE1503 AFBTSOM9358
AFBTND222AFBTND244
AFBTND255
AFBTLAG9441
AFBTND221
AFBTLAG9415AFBTND215
Eigenvalues AFBTLAG9422AFBTSOM9385
AFBTLAG9452
AFBTLAG9418
AFBTLAG9449
AFBTND211
AFBTLAG9420 AFBTLAG9419
AFBTLAG9440
AFBTLAG9445AFBTLAG9406
AFBTLAG9403
AFBTLAG9438
AFBTLAG9447
AFBTLAG9425
AFBTLAG9404
AFBTLAG9429
AFBTLAG9411AFBTLAG9428
AFBTLAG9417
AFBTLAG9450
AFBTLAG9448
AFBTLAG9421
AFBTLAG9444
AFBTLAG9426
AFBTLAG9430
AFBTLAG9431
AFBTLAG9427
AFBTLAG9435
AFBTLAG9423AFBTLAG9407
AFBTLAG9414
AFBTLAG9413
AFBTLAG9410
AFBTLAG9439
AFBTLAG9408
AFBTLAG9416
AFBTLAG9402
AFBTLAG9434
AFBTLAG9424 AFBTLAG9446
AFBTLAG9443
AFBTLAG9432
AFBTLAG9405
AFBTLAG9433
AFBTLAG9451
AFBTLAG9436
AFBTLAG9437
AFBTLAG9442
AFBTLAG9412

53
However, this figure can largely be improved. First, we can use [Link] to represent both
the genotypes and inertia ellipses for populations.

[Link](pca1$li, pop(microbov))
title("PCA of microbov dataset\naxes 1-2")
[Link](pca1$eig[1:20], 3,1,2)

d=1
PCA of microbov dataset
axes 1−2

Zebu

Borgou

Aubrac
Salers
Gascon
BlondeAquitaine
Charolais
Limousin
BretPieNoire
MaineAnjou
Montbeliard
Bazadais Somba
NDama

Eigenvalues

Lagunaire

This plane shows that the main structuring is between African an French breeds, the second
structure reflecting genetic diversity among African breeds. The third axis reflects the
diversity among French breeds:

[Link](pca1$li,pop(microbov),xax=1,yax=3,sub="PCA 1-3",csub=2)
title("PCA of microbov dataset\naxes 1-3")
[Link](pca1$eig[1:20],nf=3,xax=1,yax=3)

54
 d=1
PCA of microbov dataset
axes 1−3

Bazadais

Gascon
Aubrac
Limousin Somba
BlondeAquitaine
Salers Borgou
NDama
Montbeliard
ZebuLagunaire

BretPieNoire
Charolais

Eigenvalues
MaineAnjou

PCA 1−3
Overall, all breeds seem well differentiated.

However, we can yet improve these scatterplots, which are fortunately easy to customize.
For instance, we can remove the grid, choose different colors for the groups, use larger dots
and transparency to better assess the density of points, and remove internal segments of the
ellipses:

col <- funky(15)

[Link](pca1$li, pop(microbov),xax=1,yax=3, col=transp(col,.6), axesell=FALSE,
cstar=0, cpoint=3, grid=FALSE)

55
 Bazadais

Gascon
Aubrac
Limousin Somba
BlondeAquitaine
Salers Borgou
NDama
Montbeliard
ZebuLagunaire

BretPieNoire
Charolais

MaineAnjou

Let us now assume that we ignore the group memberships. We can still use color in an
informative way. For instance, we can recode the principal components represented in the
scatterplot on the RGB scale:

colorplot(pca1$li, pca1$li, transp=TRUE, cex=3, xlab="PC 1", ylab="PC 2")

title("PCA of microbov dataset\naxes 1-2")
abline(v=0,h=0,col="grey", lty=2)

56
 PCA of microbov dataset
axes 1−2

2
1
PC 2

0
−1
−2

−2 −1 0 1 2

PC 1

Colors are based on the first three PCs of the PCA, recoded respectively on the red, green,
and blue channel. In this figure, the genetic diversity is represented in two complementary
ways: by the distances (further away = more genetically different), and by the colors (more
different colors = more genetically different).

We can represent the diversity on the third axis similarly:

colorplot(pca1$li[c(1,3)], pca1$li, transp=TRUE, cex=3, xlab="PC 1", ylab="PC 3")

title("PCA of microbov dataset\naxes 1-3")
abline(v=0,h=0,col="grey", lty=2)

57
 PCA of microbov dataset
axes 1−3

2
1
0
PC 3

−1
−2

−2 −1 0 1 2

PC 1

6.3 Performing a Correspondance Analysis on genpop objects

Being contingency tables, the @tab slot in genpop objects can be submitted to a
Correspondance Analysis (CA) to seek a typology of populations. The approach is very
similar to the previous one for PCA.

data(microbov)
obj <- genind2genpop(microbov)

##
## Converting data from a genind to a genpop object...
##
## ...done.

ca1 <- [Link](tab(obj),scannf=FALSE,nf=3)

barplot(ca1$eig,main="Correspondance Analysis eigenvalues",
col=[Link](length(ca1$eig)))

58
 Correspondance Analysis eigenvalues

0.25
0.20
0.15
0.10
0.05
0.00

Now we display the resulting typology using a basic scatterplot:

[Link](ca1$li, sub="CA 1-2",csub=2)

[Link](ca1$eig,nf=3,xax=1,yax=2,posi="bottomright")

59
 d = 0.5

Zebu

Borgou
Salers
Aubrac

BlondeAquitaine
MaineAnjou
Gascon
Limousin
Charolais
Montbeliard
BretPieNoire
Bazadais

NDama
Somba
Eigenvalues

Lagunaire

CA 1−2
The same graph is derived for the axes 2-3:

[Link](ca1$li,xax=2,yax=3,lab=popNames(obj),sub="CA 1-3",csub=2)
[Link](ca1$eig,nf=3,xax=2,yax=3,posi="topleft")

60
 Eigenvalues d = 0.5

MaineAnjou

Charolais
BretPieNoire

Zebu
Lagunaire NDama Borgou
Somba Montbeliard
BlondeAquitaine
Gascon
Salers
Limousin
Aubrac
Bazadais

CA 1−3
As in the PCA above, axes are to be interpreted separately in terms of continental
differentiation, and between-breeds diversity. Importantly, as in any analysis carried out
at a population level, all information about the diversity within populations is lost in this
analysis. See the tutorial on DAPC for an individual-based approach which is nontheless
optimal in terms of group separation (dapc).

Note that as an alternative, wordcloud can be used to avoid overlaps in labels:

library(wordcloud)
## Error in library(wordcloud): there is no package called ’wordcloud’
[Link](1)
[Link](ca1$li*1.2, sub="CA 1-2",csub=2, clab=0, cpoint="")
textplot(ca1$li[,1], ca1$li[,2], words=popNames(obj),
cex=1.4, new=FALSE, xpd=TRUE)
## Error in textplot(ca1$li[, 1], ca1$li[, 2], words = popNames(obj), cex =
1.4, : could not find function "textplot"
[Link](ca1$eig,nf=3,xax=1,yax=2,posi="bottomright")

61
 d = 0.5

Eigenvalues

CA 1−2
However, only general trends can be interpreted: labels positions are randomised to avoid
overlap, so they no longer accurately position populations on the factorial axes.

62
7 Spatial analysis
The R software probably offers the largest collection of spatial methods among statistical
software. Here, we briefly illustrate two methods commonly used in population genetics.
Spatial multivariate analysis is covered in a dedicated tutorial; see spca tutorial for more
information.

7.1 Isolation by distance

7.1.1 Testing isolation by distance
Isolation by distance (IBD) is tested using Mantel test between a matrix of genetic distances
and a matrix of geographic distances. It can be tested using individuals as well as populations.
This example uses cat colonies from the city of Nancy. We test the correlation between
Edwards’ distances and Euclidean geographic distances between colonies.

data(nancycats)
toto <- genind2genpop(nancycats)

##
## Converting data from a genind to a genpop object...
##
## ...done.

Dgen <- [Link](toto,method=2)

Dgeo <- dist(nancycats$other$xy)
ibd <- [Link](Dgen,Dgeo)
ibd

## Monte-Carlo test
## Call: [Link](m1 = Dgen, m2 = Dgeo)
##
## Observation: -0.07011224
##
## Based on 999 replicates
## Simulated p-value: 0.722
## Alternative hypothesis: greater
##
## [Link] Expectation Variance
## -0.640766688 -0.002400549 0.011166768

plot(ibd)

63
 Histogram of sim

150
100
Frequency

50
0

−0.4 −0.2 0.0 0.2 0.4

sim

The original value of the correlation between the distance matrices is represented by the dot,
while histograms represent permuted values (i.e., under the absence of spatial structure).
Significant spatial structure would therefore result in the original value being out of the
reference distribution. Here, isolation by distance is clearly not significant.

Let us provide another example using a dataset of individuals simulated under an IBD
model:

data(spcaIllus)
x <- spcaIllus$dat2B
Dgen <- dist(x$tab)
Dgeo <- dist(other(x)$xy)
ibd <- [Link](Dgen,Dgeo)
ibd

## Monte-Carlo test
## Call: [Link](m1 = Dgen, m2 = Dgeo)
##
## Observation: 0.1267341

64
##
## Based on 999 replicates
## Simulated p-value: 0.001
## Alternative hypothesis: greater
##
## [Link] Expectation Variance
## 3.493743e+00 3.279537e-05 1.315165e-03

plot(ibd)

Histogram of sim
200
150
Frequency

100
50
0

−0.10 −0.05 0.00 0.05 0.10 0.15

sim

This time there is a clear isolation by distance pattern.

7.1.2 Cline or distant patches?

The correlation between genetic and geographic distances can occur under a range of different
biological scenarios. Classical IBD would result in continuous clines of genetic differentiation
and cause such correlation. However, distant and differentiated populations would also result
in such a pattern. These are slightly different processes and we would like to be able to
disentangle them. A very simple first approach is simply plotting both distances:

65
plot(Dgeo, Dgen)
dist_lm <- lm([Link](Dgen) ~ [Link](Dgeo))
abline(dist_lm, col="red", lty=2)

7
6
Dgen

5
4

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

Dgeo

Most of the time, simple scatterplots fail to provide a good picture of the data as the density
of points in the scatterplot is badly displayed. Colors can be used to provide better (and
prettier) plots. Local density is measured using a 2-dimensional kernel density estimation
(kde2d), and the results are displayed using image; colorRampPalette is used to generate a
customized color palette:

library(MASS)
dens <- kde2d([Link](Dgeo), [Link](Dgen), n=300)
myPal <- colorRampPalette(c("white", "blue", "gold", "orange", "red"))
plot(Dgeo, Dgen, pch=20,cex=.5)
image(dens, col=transp(myPal(300),.7), add=TRUE)
abline(dist_lm)
title("Isolation by distance plot")

66
 Isolation by distance plot

7
6
Dgen

5
4

0.0 0.2 0.4 0.6 0.8 1.0 1.2 1.4

Dgeo

The scatterplot clearly shows one single consistent cloud of point, without discontinuities
which would have indicated patches. This is reassuring, since the data were actually simulated
under an IBD (continuous) model.

7.2 Using Monmonier’s algorithm to define genetic boundaries

Monmonier’s algorithm [9] was originally designed to find boundaries of maximum differences
between contiguous polygons of a tesselation. As such, the method was basically used in
geographical analysis. More recently, [8] suggested that this algorithm could be employed to
detect genetic boundaries among georeferecend genotypes (or populations). This algorithm
is implemented using a more general approach than the initial one in adegenet.
Instead of using Voronoi tesselation as in the original version, the functions monmonier
and [Link] can handle various neighbouring graphs such as Delaunay
triangulation, Gabriel’s graph, Relative Neighbours graph, etc. These graphs define spatial
connectivity among locations (of genotypes or populations), with couple of locations
being neighbours (if connected) or not. Another information is given by a set of markers
which define genetic distances among these ’points’. The aim of Monmonier’s algorithm
is to find the path through the strongest genetic distances between neighbours. A more
complete description of the principle of this algorithm will be found in the documentation of

67
monmonier. Indeed, the very purpose of this tutorial is simply to show how it can be used
on genetic data.

Let’s take the example from the function’s manpage and detail it. The dataset used is
sim2pop.

data(sim2pop)
sim2pop

## /// GENIND OBJECT /////////

##
## // 130 individuals; 20 loci; 241 alleles; size: 192.1 Kb
##
## // Basic content
## @tab: 130 x 241 matrix of allele counts
## @[Link]: number of alleles per locus (range: 7-17)
## @[Link]: locus factor for the 241 columns of @tab
## @[Link]: list of allele names for each locus
## @ploidy: ploidy of each individual (range: 2-2)
## @type: codom
## @call: old2new(object = sim2pop)
##
## // Optional content
## @pop: population of each individual (group size range: 30-100)
## @other: a list containing: xy

summary(sim2pop$pop)

## P01 P02
## 100 30

temp <- sim2pop$pop

levels(temp) <- c(3,5)
temp <- [Link]([Link](temp))
plot(sim2pop$other$xy,pch=temp,cex=1.5,xlab='x',ylab='y')
legend("topright",leg=c("Pop A", "Pop B"),pch=c(3,5))

68
 100
Pop A
Pop B
80
60
y

40
20
0

0 20 40 60 80 100

There are two sampled populations in this dataset, with inequal sample sizes (100 and 30).
Twenty microsatellite-like loci are available for all genotypes (no missing data). monmonier
requires several arguments to be specified:

args(monmonier)

## function (xy, dist, cn, threshold = NULL, [Link] = NULL, nrun = 1,

## [Link] = rep(0, nrun), scanthres = [Link](threshold),
## allowLoop = TRUE)
## NULL
The first argument (xy) is a matrix of geographic coordinates, already stored in sim2pop.
Next argument is an object of class dist, which is the matrix of pairwise genetic distances.
For now, we will use the classical Euclidean distance between allelic profiles of the individuals.
This is obtained by:

D <- dist(sim2pop$tab)

The next argument (cn) is a connection network. Routines for building such networks are
scattered over several packages, but all made available through the function chooseCN. Here,

69
we disable the interactivity of the function (ask=FALSE) and select the second type of graph
which is the graph of Gabriel (type=2).

gab <- chooseCN(sim2pop$other$xy, ask=FALSE, type=2)

## Registered S3 method overwritten by ’spdep’:

## method from
## [Link] ape

The obtained network is automatically plotted by the function. It seems we are now ready
to proceed to the algorithm.

mon1 <- monmonier(sim2pop$other$xy, D, gab)

70
This plot shows all local differences sorted in decreasing order. The idea behind this is that
a significant boundary would cause local differences to decrease abruptly after the boundary.
This should be used to choose the threshold difference for the algorithm to stop extending
the boundary. Here, there is no indication af an actual boundary.

Why do the algorithm fail to find a boundary? Either because there is no genetic
differentiation to be found, or because the signal differentiating both populations is too
weak to overcome the random noise in genetic distances. What is the Fst between the two
samples?

library(hierfstat)
[Link](sim2pop, method = "Nei87")

## P01
## P02 0.0685

This value would be considered as very weak differentiation (FST = 0.023). Is it significant?
We can easily ellaborate a permutation test of this FST value; to save computational time,
we use only a small number of replicates to generate FST values in absence of population
structure:

replicate(10, {
pop(sim2pop) <- sample(pop(sim2pop))
[Link](sim2pop, method = "Nei87")
})

## [1] -0.0001 0.0014 0.0014 0.0007 -0.0023 0.0004 -0.0037 -0.0027 -0.0005 -0.0031

71
 FST values in absence of population structure would be one order of magnitude lower
(more replicate would give a very low p-value — just replace 10 by 200 in the above
command). In fact, the two samples are indeed genetically differentiated.

Can Monmonier’s algorithm find a boundary between the two populations? Yes, if we
get rid of the random noise. This can be achieved using a simple ordination method such as
Principal Coordinates Analysis.

pco1 <- [Link](D, scannf=FALSE, nf=1)

barplot(pco1$eig, main="Eigenvalues")

Eigenvalues
2.0
1.5
1.0
0.5
0.0

We retain only the first eigenvalue. The corresponding coordinates are used to redefine the
genetic distances among genotypes. The algorithm is then re-run.

D <- dist(pco1$li)

72
mon1 <- monmonier(sim2pop$other$xy, D, gab)

mon1

##
## ###########################################################
## # List of paths of maximum differences between neighbours #
## # Using a Monmonier based algorithm #
## ###########################################################
##
## $call:monmonier(xy = sim2pop$other$xy, dist = D, cn = gab, scanthres = FALSE)
##
## # Object content #
## Class: monmonier
## $nrun (number of successive runs): 1
## $run1: run of the algorithm
## $threshold (minimum difference between neighbours): 1.630755
## $xy: spatial coordinates
## $cn: connection network
##
## # Runs content #
## # Run 1
## # First direction
## Class: list
## $path:
## x y

73
## Point_1 14.98299 93.81162
##
## $values:
## 4.563555
## # Second direction
## Class: list
## $path:
## x y
## Point_1 14.98299 93.81162
## Point_2 30.74508 87.57724
## Point_3 33.66093 86.14115
## ...
##
## $values:
## 4.563555 3.23581 3.906439 ...

This may take some time... but never more than five minutes on an ’ordinary’ personnal
computer. The object mon1 contains the whole information about the boundaries found. As
several boundaries can be seeked at the same time (argument nrun), you have to specify
about which run and which direction you want to get informations (values of differences or
path coordinates). For instance:

names(mon1)

## [1] "run1" "nrun" "threshold" "xy" "cn" "call"

names(mon1$run1)

## [1] "dir1" "dir2"

mon1$run1$dir1

## $path
## x y
## Point_1 14.98299 93.81162
##
## $values
## [1] 4.563555

It can also be useful to identify which points are crossed by the barrier; this can be done
using [Link]:

[Link](mon1)

## $run1

74
## $run1$dir1
## [Link] [Link] first second
## Point_1 14.98299 93.81162 11 125
##
## $run1$dir2
## [Link] [Link] first second
## Point_1 14.98299 93.81162 11 125
## Point_2 30.74508 87.57724 44 128
## Point_3 33.66093 86.14115 20 128
## Point_4 35.28914 81.12578 68 128
## Point_5 33.85756 74.45492 68 117
## Point_6 38.07622 71.47532 68 122
## Point_7 41.97494 70.02783 35 122
## Point_8 43.45812 67.12026 69 122
## Point_9 42.20206 59.59613 22 122
## Point_10 42.48613 52.55145 22 124
## Point_11 40.08702 48.61795 13 124
## Point_12 39.20791 43.89978 13 127
## Point_13 38.81236 40.34516 62 127
## Point_14 37.32112 36.35265 62 130
## Point_15 37.96426 30.82105 94 130
## Point_16 32.79703 28.00517 16 130
## Point_17 30.12832 28.60376 85 130
## Point_18 20.92496 29.21211 63 119
## Point_19 16.05811 22.72600 61 126
## Point_20 11.72524 21.15519 89 126
## Point_21 10.18696 16.61536 74 89

The returned dataframe contains, in this order, the x and y coordinates of the points of the
barrier, and the identifiers of the two ’parent’ points, that is, the points whose barycenter is
the point of the barrier.
Finally, you can plot very simply the obtained boundary using the method plot:

plot(mon1)

75
see arguments in ?[Link] to customize this representation. Last, we can compare
the infered boundary with the actual distribution of populations:

plot(mon1, [Link]=FALSE, bwd=10, col="black")

points(sim2pop$other$xy, cex=2, pch=20,
col=fac2col(pop(sim2pop), [Link]=spectral))
legend("topright", leg=c("Pop A", "Pop B"), pch=c(20),
col=spectral(2), [Link]=2)

76
 Pop A
Pop B

Not too bad...

8 Simulating hybridization
The function hybridize allows to simulate hybridization between individuals from two
distinct genetic pools, or more broadly between two genind objects. Here, we use the
example from the manpage of the function, to go a little further. Please have a look at
the documentation, especially at the different possible outputs (outputs for the software
STRUCTURE is notably available).

temp <- seppop(microbov)

names(temp)

## [1] "Borgou" "Zebu" "Lagunaire" "NDama" "Somba"

## [10] "Charolais" "Gascon" "Limousin" "MaineAnjou" "Montbel

salers <- temp$Salers

zebu <- temp$Zebu
zebler <- hybridize(salers, zebu, n=40, pop="zebler")

77
A first generation (F1) of hybrids ’zebler’ is obtained. Is it possible to perform a backcross,
say, with ’salers’ population? Yes, here it is:

F2 <- hybridize(salers, zebler, n=40)

F3 <- hybridize(salers, F2, n=40)
F4 <- hybridize(salers, F3, n=40)

Finally, note that despite this example shows hybridization between diploid organisms,
hybridize is not restrained to this case. In fact, organisms with any even level of ploidy can
be used, in which case half of the genes is taken from each reference population. Ultimately,
more complex mating schemes could be implemented... suggestion or (better) contributions
are welcome!

78
References
[1] D. Charif and J.R. Lobry. SeqinR 1.0-2: a contributed package to the R project
for statistical computing devoted to biological sequences retrieval and analysis. In
H.E. Roman U. Bastolla, M. Porto and M. Vendruscolo, editors, Structural approaches
to sequence evolution: Molecules, networks, populations, Biological and Medical Physics,
Biomedical Engineering, pages 207–232. Springer Verlag, New York, 2007. ISBN : 978-
3-540-35305-8.

[2] S. Dray and A.-B. Dufour. The ade4 package: implementing the duality diagram for
ecologists. Journal of Statistical Software, 22(4):1–20, 2007.

[3] T. Jombart. adegenet: a R package for the multivariate analysis of genetic markers.
Bioinformatics, 24:1403–1405, 2008.

[4] T. Jombart, S. Devillard, A.-B. Dufour, and D. Pontier. Revealing cryptic spatial
patterns in genetic variability by a new multivariate method. Heredity, 101:92–103,
2008.

[5] T. Jombart, R. M. Eggo, P. J. Dodd, and F. Balloux. Reconstructing disease outbreaks

from genetic data: a graph approach. Heredity, 106:383–390, 2010.

[6] T. Jombart, D. Pontier, and A-B. Dufour. Genetic markers in the playground of
multivariate analysis. Heredity, 102:330–341, 2009.

[7] Thibaut Jombart, Sebastien Devillard, and Francois Balloux. Discriminant analysis
of principal components: a new method for the analysis of genetically structured
populations. BMC Genetics, 11(1):94, 2010.

[8] F. Manni, E. Guérard, and E. Heyer. Geographic patterns of (genetic, morphologic,

linguistic) variation: how barriers can be detected by ”Monmonier’s algorithm”. Human
Biology, 76:173–190, 2004.

[9] M. Monmonier. Maximum-difference barriers: an alternative numerical regionalization

method. Geographical Analysis, 3:245–261, 1973.

[10] M. Nei. Analysis of gene diversity in subdivided populations. Proc Natl Acad Sci U S
A, 70(12):3321–3323, Dec 1973.

[11] R Development Core Team. R: A Language and Environment for Statistical Computing.
R Foundation for Statistical Computing, Vienna, Austria, 2011. ISBN 3-900051-07-0.

[12] B. S. Weir and C. C. Cockerham. Estimating f -statistics for the analysis of population
structure. Evolution, 38:1350–1370, 1984.

79