4 - Landscape and genetic graph comparison with `graph4lg`"

knitr::opts_chunk$set(
collapse = TRUE,
comment = "#>"
)

library(graph4lg)
library(igraph)

Introduction

The rationale of graph4lg package in R is to make easier the construction and analysis of genetic and landscape graphs in landscape genetic studies (hence the name graph4lg, meaning Graphs for Landscape Genetics). This package provides users with tools for:

Each one of the included tutorials focuses on one of these points. This fourth and last tutorial will then focus on landscape and genetic graph comparison. It will describe the package functions allowing users to:

Data used in this tutorial

The package already includes genetic and spatial simulated data sets allowing users to discover its different functionalities. The first data set (data_simul) was simulated 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. 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.

The second data set (data_ex) was simulated as the first one but included only 10 populations. It is used to generate quick examples.

Here, we also rely on a data set created only for the vignettes (data_tuto) and containing several objects:

data("data_tuto")

mat_dps <- data_tuto[[1]]
mat_pg <- data_tuto[[2]]
graph_ci <- data_tuto[[3]]
dmc <- data_tuto[[4]]
land_graph <- data_tuto[[5]]
mat_ld <- data_tuto[[6]]

We will also use Graphab projects already created for the vignettes.

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, we can compare landscape and genetic graphs at these different scales, in particular when they share the same nodes. We present how to implement these comparisons using graph4lg. To see how to create landscape and genetic graphs with graph4lg, we invite users to read the second and third tutorials included in this package.

Comparing node metrics

First, landscape and genetic graphs can be compared by comparing connectivity metrics measured at the level of a habitat patch (landscape graph node) with the genetic response of the population living and sampled in this habitat patch (genetic graph node) in terms of genetic diversity and differentiation from the other populations. When using this approach, two graphs are similar if the correlation coefficient between the metric values calculated for the same nodes is high. This correlation would be interpreted as an evidence for the influence of habitat connectivity on population genetic structure.

The function graph_node_compar computes these correlations to compare two graphs at the node level. It takes as arguments:

In order to include node attributes in the correlation, users can use compute them and associate them with graph nodes with add_nodes_attr function (described in tutorial 2).

In the following example, we will import a landscape graph made of the habitat patches occupied by the 50 populations used in a gene flow simulation with CDPOP (pts_pop_simul and data_simul_genind data set). A genetic graph will be created from the simulated genetic data set. We will then compare two metrics computed at the node-level in these graphs.

land_graph is the landscape graph made of the forest habitat patches occupied by the 50 populations. Its 1225 links are weighted by cost-distances between patches (mat_ld). It is a complete graph.

land_graph <- gen_graph_topo(mat_w = mat_ld,
                             mat_topo = mat_ld,
                             topo = "comp")

# Plot the histogram of its link weights
plot_w_hist(graph = land_graph)

Using the function plot_w_hist, we see that its link weights almost follow a normal distribution with values ranging from 0 to 10.000 cost-distance units.

We will compute the mean of the inverse weight of every link connected to its nodes. This metric miw is akin to the Flux metric computed in Graphab. It takes high values if a patch is well connected and close to other patches. We use the function compute_node_metric

miw_lg <- compute_node_metric(graph = land_graph, metrics = "miw")
head(miw_lg)

We include the values of this metric as a graph node attribute with the function add_nodes_attr:

land_graph <- add_nodes_attr(graph = land_graph,
                             input = "df",
                             data = miw_lg,
                             index = "ID")

We now create a genetic graph from the simulated data set. To that purpose, we first compute a genetic distance matrix using the D~PS~ genetic distance.

mat_dps <- mat_gen_dist(x = data_simul_genind, dist = "DPS")

We use this distance matrix to create a complete genetic graph with gen_graph_topo:

gen_comp_graph <- gen_graph_topo(mat_w = mat_dps,
                                 mat_topo = mat_dps,
                                 topo = "comp")

We plot the distribution of the link weights.

plot_w_hist(graph = gen_comp_graph, 
            fill = "darkblue")

We will compute the mean of the inverse weights of links connected to each node in the genetic graph gen_comp_graph and include it as a node attribute.

miw_comp <- compute_node_metric(graph = gen_comp_graph, metrics = "miw")
gen_comp_graph <- add_nodes_attr(graph = gen_comp_graph,
                                 input = "df",
                                 data = miw_comp,
                                 index = "ID")

We can now assess the correlation between these two metrics using graph_node_compar:

graph_node_compar(x = land_graph, y = gen_comp_graph,
                  metrics = c("miw", "miw"), method = "spearman",
                  weight = TRUE, test = TRUE)

The two metrics are highly correlated, meaning that if a habitat patch is well connected, the population occupying this patch is less different from the other populations from a genetic point of view than if the habitat patch is not connected well. This is quite an expected result given that gene flow was largely driven by cost-distances between populations during the genetic simulations.

Using complete graphs, visual representations are unclear because of link overlap. Building a Gabriel graph is a way to simplify graph topology. Gabriel graphs links will be computed from the Euclidean geographical distances between populations.

mat_geo <- mat_geo_dist(data = pts_pop_simul,
                        ID = "ID", x = "x", y = "y")
mat_geo <- reorder_mat(mat_geo, order = row.names(mat_dps))

gen_gab_graph <- gen_graph_topo(mat_w = mat_dps,
                                mat_topo = mat_geo,
                                topo = "gabriel")

# Associate the values of miw from the complete graph to this graph

gen_gab_graph <- add_nodes_attr(gen_gab_graph,
                                data = miw_comp,
                                index = "ID")

# Plot the graph with node sizes proportional to MIW

plot_graph_lg(graph = gen_gab_graph, 
              crds = pts_pop_simul,
              mode = "spatial",
              node_size = "miw",
              link_width = "inv_w")

Comparing graph links

Another way to compare genetic and landscape graphs is to compare their links, thereby focusing on their topology. Two graphs are similar if the conserved links of one graph are the same ones as those conserved in the other. 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, we 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.

The function takes as arguments:

For example, we will create a landscape graph from the 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.

We first create the thresholded landscape graph and plot it

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

plot_graph_lg(land_graph_thr, 
              mode = "spatial", 
              crds = pts_pop_simul,
              link_width = "inv_w", 
              pts_col = "#80C342")

We then compare the topology of land_graph_thr and gen_gab_graph using the Matthews' correlation coefficient using graph_topo_compar:

graph_topo_compar(obs_graph = land_graph_thr, 
                  pred_graph = gen_gab_graph,
                  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, we 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. It takes as arguments:

For example:

graph_plot_compar(x = land_graph_thr, y = gen_gab_graph,
                  crds = pts_pop_simul)

This representation clearly indicates where are the shared and unshared links.

Comparing graph modules

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.

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 node pairs classified in the same module in both graphs. The function also performs the nodes partition into modules by using the modularity calculations available in igraph package. We can specify the algorithm used to compute the modularity, the link weighting, the number of modules, among others. It takes as arguments:

Two different weightings can be used to create the modules of the two graphs:

In the following example, we compare land_graph_thr and gen_gab_graph with the default parameters (algo='fast_greedy' algorithm, optimal number of modules, links weighted by inverse distances (node_inter="distance")).

graph_modul_compar(x = land_graph_thr, 
                   y = gen_gab_graph)

The ARI value is relatively high meaning that modules are somehow similar.

We can check this result visually. We first compute the modules in both graphs and include them as node attributes.

module_land <- compute_graph_modul(graph = land_graph_thr, 
                                   algo = "fast_greedy",
                                   node_inter = "distance")

land_graph_thr <- add_nodes_attr(graph = land_graph_thr,
                                 data = module_land,
                                 index = "ID")

module_gen <- compute_graph_modul(graph = gen_gab_graph, 
                                   algo = "fast_greedy",
                                   node_inter = "distance")

gen_gab_graph <- add_nodes_attr(graph = gen_gab_graph,
                                 data = module_gen,
                                 index = "ID")

We then plot both graphs, colouring their nodes depending on their module:

plot_graph_lg(graph = land_graph_thr,
              mode = "spatial",
              crds = pts_pop_simul,
              module = "module")
plot_graph_lg(graph = gen_gab_graph,
              mode = "spatial",
              crds = pts_pop_simul,
              module = "module")

We easily visualise divergence and convergence in both module partitions.

Conclusion

graph4lg makes available a large range of tools to perform landscape genetic analyses relying upon genetic and landscape graphs. We hope that next improvements will draw on results obtained using these tools in empirical as well as simulation studies.

References



Try the graph4lg package in your browser

Any scripts or data that you put into this service are public.

graph4lg documentation built on March 16, 2021, 5:09 p.m.