---
title: '`graph4lg` - Tutorial'
author: "Paul SAVARY"
date: "`r Sys.Date()`"
output:
  html_document:
    df_print: paged
bibliography: biblio_vignette.bib
vignette: |
  %\VignetteIndexEntry{`graph4lg` - Tutorial} 
  %\VignetteEngine{knitr::rmarkdown} 
  %\VignetteEncoding{UTF-8}
---
  
```{r setup, include = FALSE}
knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)

library(graph4lg)
library(igraph)

```

## Introduction

The creation of the `graph4lg` package on R should make easier the construction and analysis of genetic and landscape graphs in landscape genetics studies (hence the name `graph4lg`). The package allows to weight the links and to prune the graphs by several ways. To our knowledge, it is the first software which enables to create genetic graphs with such a large variety of parameters. Besides, it allows to carry out preliminary analyses of the spatial pattern of genetic differentiation in order to choose the best genetic distance and pruning method in every context. Lastly, it makes possible the comparison of two spatial graphs sharing the same nodes.


In the following sections, we describe the functions of the package following their separation into three parts. We illustrate them with results' examples and report the command lines to copy in order to reproduce them.

## Loading of the data used to create the tutorial

The package includes genetic and spatial data allowing users to discover its different functionalities. We included two datasets. The first is from French Guadeloupe and consists of 432 Plumbeous Warblers (*Setophaga plumbea*) genotyped at 12 microsatellite loci (169 different alleles in total) from 20 populations whose spatial coordinates are available. These data were collected and analysed by @khimoun2017landscape. They are also freely available on the website [datadryad.org](http://datadryad.org/). They were included in the package in different formats: `data.frame` format with columns of `locus` class from the package `gstud`, `genind` class from `adegenet`, `loci` class for `pegas`, text files (.txt) for GENEPOP or STRUCTURE.

The second data set was created during simulations done with CDPOP [@landguth2010cdpop] on a simulated landscape. It consists of 1500 individuals from 50 populations genotyped at 20 microsatellite loci. Individuals dispersed less when the cost-distance between populations was large. These data were included in format `genind` and as a text file compatible with GENEPOP software. A landscape graph was created with GRAPHAB [@foltete2012software] whose nodes were the 50 simulated populations and the links were weighted by cost-distance values between populations. The project created with GRAPHAB was included into the package such that the landscape graphs and the cost-distance matrix can be easily imported into the R environment.


```{r, echo = FALSE, eval = TRUE}
data("data_tuto")

mat_dps <- data_tuto[[1]]
mat_pg <- data_tuto[[2]]
graph_ci <- data_tuto[[3]]
dmc <- data_tuto[[4]]

```


## Input data processing

The first type of functions from this package allows to process spatial and genetic data used as input data in subsequent analyses performed with other functions. These functions convert genetic data, compute genetic or geographical distances, import landscape graphs and convert cost-distance values into Euclidean geographical distances.

### Genetic data conversion


In order to make the package user-friendly and compatible with genetic data commonly used in landscape genetics, the functions `gstud_to_genind`, `loci_to_genind`, `structure_to_genind` and `genepop_to_genind` allow to convert genetic data from formats used respectively in the GENETIC STUDIO [@dyer2009geneticstudio] and `pegas` [@paradis2010pegas] packages on R and in STRUCTURE [@pritchard2000inference] and GENEPOP [@raymond1995genepop] software into R objects with the class attribute `genind` from ADEGENET package [@jombart2008adegenet]. The format `genind` allows to use time-efficient functions from ADEGENET (coded in C). This package was developed and is regularly maintained by Thibaut Jombart [(his website)](https://thibautjombart.netlify.com/). The function `genind_to_genepop` enables to convert `genind` object into text files in format `genepop` in order to perform easily analyses with this commonly used R package and executable software.

Genetic data in formats used by software commonly used in population genetics can be converted into `genind` objects in R. This class of object was created for the package `adegenet` [@jombart2008adegenet]. 



#### GENEPOP to genind

The GENEPOP software [@raymond1995genepop] developed by M. Raymond and F. Rousset (Montpellier) can be used as an executable file, with or without graphical user interface, or as a R package. It is frequently used to compute F~ST~ values and to test for Hardy-Weinberg equilibrium, linkage disequilibrium or genetic differentiation.


When performing simulations with CDPOP, individuals' genotypes can be saved as GENEPOP files at the end of the simulation.

The function `genepop_to_genind` loads a GENEPOP file (.txt extension) and converts it into a `genind` object. To use it, the path to the file, the total number and names of the loci and the populations' names must be indicated.

```{r}
data_genind <- genepop_to_genind(path = paste0(system.file('extdata', 
                                                           package = 'graph4lg'), "/gpop_51_sim22_01_25.txt"),
                                 n.loci = 20, pop_names = as.character(1:50))
data_genind
```

We get a `genind` object. It contains the genotypes of the 1500 individuals from the 50 populations created during the simulation.


#### genind to GENEPOP

Similarly, the function `genind_to_genepop` allows to perform the reverse conversion, i.e. to convert a `genind` object into a GENEPOP file. This file is created and saved in the working directory defined earlier.


```{r, echo = TRUE, eval = FALSE}
genind_to_genepop(x = data_genind, output = "data_gpop_test.txt")
```

This function allows for example to create a GENEPOP file to test for between populations genetic differentiation or to compute fixation indices with GENEPOP software.

#### STRUCTURE to genind

STRUCTURE software [@pritchard2000inference] is frequently used in population genetics and landscape genetics. It enables to create populations' clusters via a Bayesian approach aiming at minimizing the deviation from Hardy-Weinberg equilibrium when gathering populations with one another. The input files have a particular structure. The function `structure_to_genind` allows to convert this type of file into a `genind` object.

To use the function, we need to indicate the path to the file, the names of the loci, the individuals' ID and the populations' names in the same order as in the original file.

```{r}
loci_names <- c("DkiD104", "DkiD124", "DkiD102", "CAM19",
                "DkiC118", "DkiD128", "DkiB12",  "Lswmu7",
                "DkiD109", "Lswmu5",  "TG12_15", "DkiD12" )
data(data_pc_genind)
ind_names <- row.names(data_pc_genind@tab)
pop_names <- c("BT-1", "BT-10", "BT-11", "BT-12", "BT-13", "BT-2",
               "BT-3",  "BT-4",  "BT-5",  "BT-6",  "BT-7",  "BT-8",
               "BT-9",  "GT-1",  "GT-2", "GT-3",  "GT-4",  "GT-5",
               "GT-6", "GT-7")
data_paru <- structure_to_genind(path = paste0(system.file('extdata', 
                                                           package = 'graph4lg'), 
                                               "/data_PC_str.txt"),
                                 loci_names = loci_names,
                                 pop_names = pop_names,
                                 ind_names = ind_names)
data_paru
```

#### gstud to genind

Packages `gstudio` and `popgraph` developed by R. Dyer [@dyer2009geneticstudio] use as input data R `data.frames` with columns of class `locus`. These `data.frame` objects constitute `gstud` objects. Given these packages are often used to create genetic graphs, we created a function to convert them into the `genind` format. 

A `gstud` object generally has the following structure (Plumbeous Warbler data from Guadeloupe as an example):
```{r}
head(data_pc_gstud)
```

To convert it with the function `gstud_to_genind`, we indicate the name of the `data.frame` and of the columns with populations' names and individuals' names:
```{r}
gstud_to_genind(x = data_pc_gstud, pop_col = "Cluster",
                ind_col = "ID")
```

### Genetic and geographical distances calculation

#### Genetic distances 

From a `genind` object, the function `mat_gen_dist` calculates several types of between populations genetic distances:

+ F~ST~ [@weir1984estimating] or linearised F~ST~ [@rousset1997genetic] (options '`dist=FST`' and '`dist=FST_lin`').

+ G'~ST~ [@hedrick2005standardized] (option '`dist=GST`').

+ D~Jost~ [@jost2008gst] (option '`dist=D`').

+ D~PS~ (1 - proportion of shared alleles) [@bowcock1994high, @murphy_graph_2015] (option '`dist=DPS`').

+ Euclidean genetic distance [@excoffier1992analysis] (option '`dist=basic`').

+ Euclidean genetic distance with a weighting depending on allelic frequencies giving more weight to rare alleles [@fortuna2009networks] (option '`dist=weight`').

+ Euclidean genetic distance computed after a PCA of the matrix of allelic frequencies by population. The axes considered to compute the Euclidean distance are the non-collinear principal components (total number of alleles - number of loci) [@paschou2014maritime, @shirk2017comparison] (option '`dist=PCA`').

+ Euclidean genetic distance computed in the same way as with the function `popgraph` from `popgraph` package, i.e. after a PCA and two SVD, among other computation steps (option '`dist=PG`'). This distance is different from the conditional genetic distance (cGD) computed from a population graph by summing genetic distances along shortest paths.

To do these calculations with the function `mat_gen_dist`, you just have to indicate the name of the `genind` object which includes the genetic data of the individuals as well as the populations to which each of them belongs. The other argument of the function is the type of genetic distance to compute. Here are two examples:

```{r,eval=FALSE, echo =TRUE, message = FALSE, warning = FALSE}
mat_dps <- mat_gen_dist(x = data_genind, dist = "DPS")
```
```{r, message = FALSE, warning = FALSE}
mat_dps[1:5, 1:5]
```
```{r, eval=FALSE, echo =TRUE, message = FALSE, warning = FALSE}
mat_pg <- mat_gen_dist(x = data_genind, dist = "PG")
```
```{r, message = FALSE, warning = FALSE}
mat_pg[1:5, 1:5]
```

#### Import of landscape graphs created with GRAPHAB

The function `graphab_to_igraph` allows to import landscape graphs created with \textsc{graphab} software into the R environment. It takes the path to the directory of a project created with \textsc{graphab} as a first argument. Then, the names of the shapefile layers corresponding to the nodes and links (the least-cost paths or a set of graph links) of the graphs must be indicated (for the nodes, by default, the shapefile layers "patches.shp" is imported). The links of the imported graph are weighted by cost-distance values or by geographical distance values depending on the `weight` option. The graph can be plotted on a geographical map. Besides, an object of type `data.frame` with the spatial coordinates of the nodes' centroids can be created.

As an example, the following command allows to import the project created to parameterize the gene-flow simulation on \textsc{cdpop} performed to create `data_simul_genind`:

```{r}
land_graph <- graphab_to_igraph(dir_path = system.file('extdata', 
                                                       package = 'graph4lg'), 
                                nodes = "patches", 
                                links = "liens_rast2_1_11_01_19-links",
                                weight = "cost", fig = FALSE, crds = TRUE)
```

The `land_graph` object is a list composed of an object of class `igraph` corresponding to the landscape graph, and of a table of class `data.frame` with the coordinates of the patches' centroids.

```{r}
crds_patches <- land_graph[[2]]
land_graph <- land_graph[[1]]
```

We create a matrix with the cost-distances between populations from the links of the landscape graph imported, which is a complete graph. You have to pay attention to the order of the populations' names in the distance matrices, even if error messages will appear when using the functions of this package if you use two matrices in which rows' and columns' names do not match. The function `reorder_mat` reorders the rows and columns of a symmetric matrix according to a specified order. To use it, you just need to create a character vector with the rows and columns names from the symmetric matrix in the order you want them to be ordered in the new matrix.

```{r}
mat_ld <- as_adjacency_matrix(land_graph, attr = "weight", 
                              type = "both", sparse = FALSE)
order <- row.names(mat_ld)[order(c(as.character(row.names(mat_ld))))]
mat_ld <- reorder_mat(mat = mat_ld, order = order)
```

#### Geographical distances

You can also calculate Euclidean geographical distances between populations with the function `mat_geo_dist`. It takes as an argument a `data.frame` with 3 columns:

+ ID : populations' ID (or simple points' ID when using this function in another context)

+ x : populations' or points' longitude

+ y : populations' or points' latitude

Geographical coordinates must be expressed in a projected coordinates system (metric) because Pythagoras's theorem is used to compute the distances (a warning message is displayed in every case).

Example :
```{r}
head(crds_patches)
```
```{r}
mat_geo <- mat_geo_dist(data = crds_patches, ID = "ID", x = "x", y = "y")
mat_geo <- reorder_mat(mat = mat_geo, order = order) 
```
```{r}
mat_geo[1:5, 1:5]
```

Cost-distances are expressed in cost units arbitrarily defined based on the cost values assigned to every land cover type when creating resistance surfaces. However, species dispersal distances are usually known as distances expressed in metric units. When we know that the probability that a species covers 10 km is 5 %, we can ask what is the equivalent of this distance in cost distance units according to the scenario of cost values assumed. It allows to prune a landscape graph with a given distance threshold e.g. or provides an order of magnitude of the cost-distance values. 

To that purpose, you can perform a regression of the between populations cost-distance values against the corresponding geographical distances. It estimates the relationship between both types of distance. Then, the resulting parameters estimates enable to convert a geographical distance into its cost-distance equivalent according to the cost scenario. 

The function `convert_cd` performs the linear regression or log-log linear regression between the geographical distance matrix and the cost-distance matrix, in the same way as @tournant2013evaluating and as performed by GRAPHAB software.

Below, we estimate the relationship between geographical distance and cost-distance between the populations used to perform the gene-flow simulation. We convert 10 km into cost-distance units. The function also plots the relationship between these distances.

```{r, out.width = '40%'}
convert_res <- convert_cd(mat_euc = mat_geo, mat_ld = mat_ld, 
                          to_convert = 10000, fig = TRUE, 
                          method = "log-log", pts_col = "grey")
convert_res
```

In this case, we can see that 10 km are equivalent to 1606 cost-distance units. The log-log linear model estimates the relationship between geographical distance (GD) and cost-distance (CD) such that: $log(CD)=-2.2512+1.0458 \times log(GD)$. The determination coefficient $R^2$ associated to this linear model is 0.69. The black dot represented on the figure refers to the 10 km value on the regression line characterizing the relationship between cost-distance and geographical distance.

## Graphs' construction

Some functions from `graph4lg` package construct graphs from either `genind` objects or genetic or landscape distance matrices. To choose a genetic distance and a pruning method for the genetic graphs' construction, we developed functions to perform preliminary analyses of the spatial pattern of genetic differentiation. Indeed, a genetic graph can be created in order to i) identify the direct dispersal paths between populations or to ii) select the set of population pairs to consider to infer landscape effects on dispersal. According to the use of a genetic graph and to the spatial pattern of genetic differentiation (type-I or type-IV pattern of IBD [@hutchison1999correlation, @van2015isolation]), the choice of a genetic distance and of a pruning method will not be the same.  

@van2015isolation computed the so-called distance of maximum correlation (DMC) as the distance between populations below which population pairs should be considered in order to maximize the correlation between landscape distance (geographical distance in their case, but applies similarly to cost-distance) and genetic distance. This distance threshold is computed by increasing iteratively the maximum distance between populations above which population pairs are not taken into account to compute the correlation. Thus, an increasing number of population pairs is considered in the inference. When the correlation coefficient between landscape distance and genetic distance reaches a maximum, the distance threshold considered is the DMC. When the DMC is equal to the maximum distance between populations, it means that an equilibrium established between gene flow and genetic drift at the scale of the study area. Conversely, when the DMC is lower than this maximum distance, it means that there is a "plateau" in the relationship between landscape distance and genetic distance because migration-drift equilibrium has not been reached yet at the scale considered. It can be due to recent modifications of the landscape which consistently reduced the connectivity in a previously connected context. In this case, graph pruning is needed to well infer landscape effect on dispersal. Similarly, genetic distances that do not assume this equilibrium should be used.

The function `dist_max_corr` calculates the DMC from two distance matrices. You need to specify the interval between two distance thresholds iteratively considered to select population pairs and compute the correlation coefficient. The function `scatter_dist`, on the other hand, allows to visualize the relationship between two distance matrices by making a scatter plot. The shape of this relationship can be compared to the four different types of IBD patterns described by @hutchison1999correlation in order to characterize the spatial pattern of genetic differentiation. 

Here are the command lines to execute to use these two functions:
```{r, eval = FALSE, echo = TRUE}
dmc <- dist_max_corr(mat_gd = mat_dps, mat_ld = mat_ld, 
                     interv = 500, pts_col = "black")
```

The `dmc` object is a list with 1) the DMC value, 2) a vector containing all the computed correlation coefficients, 3) a vector with all the distance thresholds tested and 4) a graphic object created with the `ggplot2` package.

```{r}
# DMC value
dmc[[1]]
# Correlation coefficients
dmc[[2]]
# Threshold distances tested
dmc[[3]]
```

The figure below represents the evolution of the correlation coefficient values when distance thresholds increase.


```{r, out.width='40%', eval = TRUE, echo = FALSE}

vec_t <- c(500, 1000, 1500, 2000, 2500, 3000, 3500, 4000, 4500, 5000, 5500, 6000,
           6500, 7000, 7500, 8000, 8500, 9000, 9500, 10000, 10230.05)

cc_val <- c(NA, 0.2986565, 0.3154498, 0.5188747, 0.7059633, 0.7559539, 
            0.7850267, 0.7947691, 0.8038470, 0.7853646, 0.7760106, 
            0.7641339, 0.7530264, 0.7462445, 0.7386713, 0.7333936, 
            0.7305631, 0.7226695, 0.7137972, 0.7110962, 0.7041702)


dat <- data.frame(vec_t = vec_t, cc_val = cc_val)
if(any(is.na(dat$cc_val))){
  dat <- dat[-which(is.na(dat$cc_val)), ]
}

plot_dmc <- ggplot2::ggplot(data = dat, ggplot2::aes(x = vec_t, y = cc_val)) +
  ggplot2::geom_point(color = "#999999", size = 1, shape = 16) +
  ggplot2::geom_line(color = "black") +
  ggplot2::labs(x = "Distance threshold",
                y = "Correlation coefficient") +
  ggplot2::theme_bw()
plot_dmc

```

The scatter plot is displayed with the function `scatter_dist`.

```{r, out.width='40%'}
scatter_dist(mat_gd = mat_dps, mat_ld = mat_ld, 
             pts_col = "black")
```

In this particular case, we notice a type-IV pattern of isolation by distance with a "plateau" in the relationship between cost-distance and genetic-distance (D~PS~). Graph pruning will be needed to select the population pairs to include in the inference of landscape effects on dispersal.

In the following section, we present the different pruning methods available.

### Pruning based on distance thresholds

To prune a graph whose links are weighted by distances, you can remove all the links associated to geographical or genetic distances larger (or lower) than a specific threshold distance. This distance can for example be equal to the maximum dispersal distance of an individual of the study species at the scale of its lifespan so that the resulting graph represents the direct dispersal paths of the species. It can also be equal to the DMC if the objective is to infer landscape effects on dispersal. 

The function `gen_graph_thr` takes as arguments a distance matrix used to weight the links of the resulting graph (`mat_w`) and a distance matrix on which the "thresholding" is based (`mat_thr`). The selected links are selected according to the values of this latter matrix. The argument `thr` is the numerical value of the threshold distance. If `mat_thr` is not specified, `mat_w` is used by default for the thresholding. Lastly, you have to specify if the links to remove take larger or lower values than the threshold value.

```{r}
graph_thr <- gen_graph_thr(mat_w = mat_dps, mat_thr = mat_geo,
                           thr = 12000, mode = "larger")
graph_thr
```

The function returns a graph in the form of an `igraph` object, which is consequently compatible with all functions from `igraph` package [@csardi2006igraph], one of the most used R package to create and analyse graphs (together with `sna` and `networks`). In the latter example, the graph has 50 nodes and 162 links when we prune it using a 12-km distance threshold. Its links are weighted with the values of the `mat_dps` matrix.

### Pruning based on topological constraints

A graph can be pruned according to a topological criterion. The function `gen_graph_topo` can use 4 different criteria. As with the previous function, topological criteria are applied by considering the distance values of the `mat_topo` matrix, but the links are weighted with the values of the `mat_w` matrix (except when `mat_topo` is not specified, cf. previous section).

**Gabriel graph**: in the created graph, two nodes are connected by a link if, when we draw a circle whose center is set at the middle of the segment linking them and whose radius is equal to half the length of this segment, there is no other node inside the circle. In mathematical terms, it means that there is a segment between $x$ and $y$ if and only if for every other point $z$, we have: $d_{xy}\leq \sqrt{d_{xz}^{2}+d_{yz}^{2}}$. We can compute such a graph from geographical distances [@gabriel1969new] (`graph_gab_geo` below) but also, less commonly, from genetic distances [@naujokaitis2013implications] (`graph_gab_gen` below). In the latter case, it is to some extent as if Pythagoras's theorem was applied to genetic distances, which can seem a bit strange even if this method has already been used by @naujokaitis2013implications.

```{r}
graph_gab_geo <- gen_graph_topo(mat_w = mat_dps, mat_topo = mat_geo,
                                topo = "gabriel")
graph_gab_geo
graph_gab_gen <- gen_graph_topo(mat_w = mat_dps, mat_topo = mat_dps,
                                topo = "gabriel")
```

**Minimum Spanning Tree (MST)**: it creates a minimum spanning tree, i.e a graph in which every node is connected by a link to at least another node and whose total links' weight is minimum. By definition, its number of links is equal to the number of nodes - 1.

```{r}
graph_mst <- gen_graph_topo(mat_w = mat_dps, mat_topo = mat_dps,
                            topo = "mst")
graph_mst
```

**"Percolation" graph**: the graph is created by removing iteratively some links, beginning with those with the highest weights until the graph breaks into more than one component. We conserve the link whose removal entails the creation of another component to obtain a connected graph. This method is also called the *edge-thinning method* [@urban2009graph]. Such a method is linked to percolation theory [@rozenfeld2008network]. The function `gen_graph_topo` indicates the number of conserved links and the weight of the link whose removal disconnects the graph (maximum link weight of the created graph).

```{r}
graph_percol <- gen_graph_topo(mat_w = mat_dps, mat_topo = mat_dps,
                               topo = "percol")
```


**Complete graph**: the function allows to create a complete graph from a distance matrix. In that case, there is no pruning and, by definition, all population pairs are connected.

```{r}
graph_comp <- gen_graph_topo(mat_w = mat_dps, mat_topo = mat_dps,
                             topo = "comp")
```

### Pruning based on the conditional independence principle

The last pruning method implemented by the `graph4lg` package is based upon the conditional independence principle. The function `gen_graph_indep` is largely inspired by the function `popgraph` created by R. Dyer [@dyer2004population], but does not need the package `popgraph` to function. Besides, as some calculations are performed with functions from the `adegenet` package (coded in C), it is faster than the original `popgraph` function. It is also more flexible than `popgraph` function given you can vary i) the way we compute genetic distances used to weight the links and to compute the covariance between populations, ii) the formula used to compute the covariance from squared distances or alternatively simple distances, iii) the statistical tolerance threshold, iv) the p-values adjustment and v) the returned objects created by the function. Without entering further into the details, here is an implementation example. 

```{r, eval = FALSE, echo = TRUE}
graph_ci <- gen_graph_indep(x = data_genind,
                            dist = "PCA",
                            cov = "sq",
                            adj = "holm")
```
```{r}
graph_ci
```

## Graphs' processing and analysis

Once the graphs have been created, you can perform calculation from them, visualize and export them. 

### Graphical application

#### Representation of the graph on a map

The function `plot_graph_lg` integrates functions from `igraph` and `ggplot2` to represent graphs on a map. Most frequently, graphs are spatial and a table with populations' coordinates must be given as an argument. It must have exactly the same structure as the table given as an argument to `mat_geo_dist` (3 columns : ID, x, y). The visual representation can make visible the links' weights by plotting the links with a width proportional to the weight (`width = "w"`) or the inverse weight (`width = "inv"`) of the links.

For example, with the graph `graph_mst` :
```{r, out.width = '60%'}
p <- plot_graph_lg(graph = graph_mst, crds = crds_patches,
                   mode = "spatial", weight = TRUE, width = "inv")
p
```

If the populations' spatial coordinates are not available, you can still display the graph on a two-dimensional plane. In that case, the nodes' positions are computed with @fruchterman1991graph algorithm to optimize the representation. With the graph `graph_mst`, we obtain:

```{r, out.width = '60%'}
p <- plot_graph_lg(graph = graph_mst, crds = crds_patches,
                   mode = "aspatial", weight = TRUE, width = "inv")
p
```

#### Representation of the graph's modules on a map

You can also plot the results of a graph's partition carried out by computing a modularity index with the function `plot_graph_modul`. This function is similar to `plot_graph_lg`. The main difference is that the graph's nodes from the same modules are displayed in the same color. 

Several algorithms can be used: `fast greedy` [@clauset2004finding], `louvain` [@blondel2008fast], `optimal` [@brandes2008modularity] and `walktrap` [@pons2006computing]. The number of created modules in each graph is adjustable but can be fixed depending on the optimal values of the modularity indices. Besides, the modularity calculation can take into account the links' weights or not.  When taken into account, the weight given to a link in the calculation will be i) proportional to or ii) inversely proportional to the genetic or landscape distance to which it is associated. 

In the following example, the graph is partitioned into 6 modules with the `fast_greedy` algorithm (default option).
```{r, out.width = '80%'}
plot_graph_modul(graph = graph_gab_geo, crds = crds_patches)
```


#### Representation of the links on a scatterplot

In landscape genetics, a graph is generally pruned from a distance matrix in which a set of distance values between population pairs or sample sites are chosen. This matrix is usually a genetic distance matrix. The relationship between these genetic distances and corresponding landscape distances (geographical or cost-distance) can be studied. When a scatterplot is created to do that (with the function `scatter_dist`), you can display the points corresponding to population pairs connected in the pruned graph in a different color. The function `scatter_dist_g` thereby allows to understand the pruning and to assess its intensity. 

In the following example, you can see that all connected population pairs from `graph_gab_geo` are separated by short landscape distances. Besides, not a single population pair is located where a plateau is reached in the relationship between genetic and landscape distances on the scatterplot.

```{r, out.width = '40%'}
scatter_dist_g(mat_y = mat_dps , mat_x = mat_ld, graph = graph_gab_geo)
```

#### Links' weights distribution

You can also create an histogram with the links' weights distribution with the function `plot_hist_w`.

```{r, out.width = '40%'}
p <- plot_w_hist(graph = graph_gab_gen)
p
```

In this case, you can vary the bars' colour of the histogram.

### Spatial layers export

Even if the function `plot_graph_lg` enables to visualize a spatial graph on a geographical plane, it is often useful to confront the populations' and links' locations to other types of spatial data. To that purpose, you can export the graph into shapefile layers in order to open them in a GIS. The graphs' nodes must have spatial coordinates. When exporting, you can choose to export only the nodes' shapefile layer, the links' shapefile layer or both. You can also compute metrics associated to the nodes (degree, betweenness, sum of inverse links weights). These metrics will be included in the attribute table of the exported nodes' shapefile layer. For the links, the attribute table contains the weights associated to every link. The function `graph_to_shp` takes also as an argument the coordinates reference system (CRS) in which the points' coordinates from the table are expressed. It will be the CRS of the created shapefile layers. The last argument is the suffix given to the shapefile layers' names beginning with "node" or "link".


```{r, eval = FALSE, echo = TRUE}
graph_to_shp(graph = graph_mst, crds = crds_patches, mode = "both",
             layer_name = "test_shp_mst", 
             dir_path = "wd",
             metrics = TRUE,
             crds_crs = "+proj=lcc +lat_1=49 +lat_2=44 +lat_0=46.5 +lon_0=3
                      +x_0=700000 +y_0=6600000 +ellps=GRS80 +units=m +no_defs")
```

Shapefile layers are created in the working directory and can be imported into a GIS.

### Metrics calculation


You can also compute metrics from graphs with functions from `igraph` package on R. For example, to compute the nodes' degree:
```{r}
igraph::degree(graph_percol)
```

You can also perform modularity analyses:
```{r}
igraph::cluster_fast_greedy(graph_thr, 
                            weights = 1/igraph::E(graph_thr)$weight)
```

Metrics' values can be associated to the nodes to which they correspond. To that purpose, 1) you can create a `data.frame` to store the values. It must have a column with the nodes' names in the same form as the nodes' names of the `igraph` object from which you compute the metric. 

```{r}
df_met <- data.frame(ID = V(graph_percol)$name)
df_met$deg <- igraph::degree(graph_percol)
df_met$modul <- igraph::cluster_fast_greedy(graph_thr, 
                                            weights = 1/igraph::E(graph_thr)$weight)$membership
```

Then, you add the metrics' values to the `igraph` object as nodes' attributes.
```{r}
graph_percol <- add_nodes_attr(graph_percol,
                               input = "df",
                               data = df_met,
                               index = "ID")
```
It is another and more flexible way to export metrics' values in the attribute table of the shapefile layers created with the function `graph_to_shp`.

You can also 2) associate metrics' values computed with GRAPHAB software to the nodes of the `igraph` object. In this case, you have to specify the path to a shapefile layer created with GRAPHAB whose attribute table contains a field with the graph's nodes names. 

```{r, eval = FALSE, echo = TRUE}
land_graph <- add_nodes_attr(land_graph,
                             input = "shp",
                             dir_path = system.file('extdata', package = 'graph4lg'),
                             layer = "patches",
                             index = "Id",
                             include = "Area")
```

## Graphs' comparison

In landscape genetics, the analysis of the link between landscape and genetic data can be performed at several scales. Indeed, *node-, link-, neighbourhood-* and *boundary-based analyses* are distinguished [@wagner2013conceptual]. Similarly, you can compare landscape and genetic graphs at several scales, in particular when they share the same nodes. Two graphs are similar if the correlation coefficient (non parametric rank correlation ) between the metric values calculated for the same nodes or between the distance values associated to the same links is high. Similarly, they are similar if the conserved links of one graph are the same ones as those conserved in the other. Finally, two graphs have similar topological and connectivity properties if their modules match, i.e. if two nodes classified together in the same module when the partition in modules is computed from one graph are also classified together in the same module when the partition is computed from the other graph. 

We included two functions allowing such comparisons. The function `graph_topo_compar` compares the topologies of two graphs sharing the same nodes. To do that, it creates a contingency table whose modalities are "Presence of a link" and "Absence of a link". We consider the topology of one graph to represent the "reality" of the dispersal flows that the other one is supposed to reproduce. In the (double entry) contingency table, when there is a 10 in the cell "Presence of a link" $\times$ "Presence of a link", it means that 10 links between the same population pairs are present in the two graphs compared. The 10 value corresponds to the number of true positives (TP). The three other values of the table are the number of false positives (FP), true negatives (TN) and false negatives (FN). From this table, you can compute several metrics often used to evaluate the performance of classification methods: the Matthews' correlation coefficient [@matthews1975comparison], the Kappa index, the False Discovery Rate, the Accuracy, the Sensitivity, the Specificity and the Precision.

For example, you can create a landscape graph from a cost-distance matrix using a threshold of 2000 cost-distance units and then compare the topology of this graph to that of the Gabriel graph created from the same nodes by computing the Matthews' correlation coefficient. 

```{r}

land_graph2 <- gen_graph_thr(mat_w = mat_ld, mat_thr = mat_ld,
                             thr = 2000, mode = "larger")

graph_topo_compar(obs_graph = land_graph2, 
                  pred_graph = graph_gab_geo,
                  mode = "mcc", 
                  directed = FALSE)
```

We get a Matthews' correlation coefficient of 0.54. This coefficient takes a value of 0 when the matches between topologies are no more frequent than by simple chance. It reaches 1 when the topologies are identical. 

Besides, you can compare the topologies of two graphs sharing the same nodes visually. To that purpose, their links are displayed on a map with a color depending on their presence in both graphs or in only one of them. The function `graph_plot_compar` can be used:

```{r, out.width = '70%'}
graph_plot_compar(x = land_graph2, y = graph_gab_geo, 
                  crds = crds_patches)

```

Studies on graph theory developed a range of metrics. Among them, some are computed at the node level. Hence, two graphs sharing the same nodes can be compared by correlating the metrics computed at the node level in both graphs. These metrics can also be correlated to other variables characterizing the nodes, such as the allelic richness. The function `graph_node_compar` allows to compute such correlations. The metrics which can be computed at the node level are: the degree, the closeness centrality, the betweenness centrality, the strength (sum of the weights of the links connected to a node), the sum or the mean of the weights of links connected to a node [@koen2016node]. Besides, other variables can be studied by associating them to the graph's nodes as attributes with the function `add_nodes_attr`. Pearson's, Kendall's and Spearman's correlation coefficients between the metrics or variables can be computed. A significance test can then be performed.  

In the following example, the betweenness centrality index is computed in both graphs. Then, Spearman's correlation coefficient is computed.

```{r}
graph_node_compar(x = graph_gab_geo, y = land_graph2,
                  metrics = c("btw", "btw"), method = "spearman",
                  weight = TRUE, test = TRUE)
```


The function `graph_modul_compar` compares the nodes partitions into modules. To do that, it computes the Adjusted Rand Index (ARI) [@hubert1985comparing], a standardised index which counts the number of nodes pairs classified in the same module in both graphs. The function also performs the nodes partition into modules by using the modularity calculations available from `igraph`. You can specify the algorithm used to compute the modularity, the link weighting, the number of modules, among others (cf. function `plot_graph_modul`).

In the following example, we compare the same graphs as previously with the default parameters ('fast_greedy' algorithm, optimal number of modules, links weighted by inverse distances).

```{r}
graph_modul_compar(x = land_graph2, y = graph_gab_geo)
```

The ARI value is relatively high. There are 6 modules in both graph partitions.


## References
