--- title: "Getting Started with landgraph" author: "Bill Peterman" output: rmarkdown::html_vignette: number_sections: true toc: true vignette: > %\VignetteIndexEntry{Getting Started with landgraph} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6, fig.height = 4.25 ) ``` # Introduction `landgraph` holds the shared inputs that landscape-genetic network methods consume: a lightweight spatial graph, genetic covariance and distance matrices built from molecular data, and antisymmetric directional edge covariates. It is the common base beneath two downstream packages: `terradish`, which estimates symmetric landscape resistance, and `dragonflow`, which estimates asymmetric (directional) gene flow. You can also use `landgraph` on its own whenever you need a graph, a genetic covariance, or a directional covariate. This vignette walks the whole pipeline end to end. You will build a graph from deme coordinates, turn molecular data into a covariance matrix and a distance matrix, read what those matrices mean, and construct the two kinds of directional edge covariate. Every output is shown and then interpreted, so you can recognize and trust what you produce on your own data. You need only a working knowledge of R and of population-genetic terms such as allele, locus, and population. Terms specific to this package are defined as they appear. The package depends on base R and `stats` only; the spatial reader (`terra`) and the Delaunay graph builder (`deldir`) are optional and used only where noted. A note on vocabulary. A **deme** is a local breeding group of individuals, the unit at a single sampled site. A **vertex** (or **node**) is a deme's position in the graph. An **edge** joins two vertices that the graph treats as neighbors. These words are used interchangeably with their graph meaning throughout. ```{r load} library(landgraph) ``` # Building a deme graph The graph is the spatial scaffold. `deme_graph()` takes a two-column matrix of coordinates, one row per deme, and returns the vertices plus an undirected edge list. We will lay six demes on a small grid and connect them with **rook** adjacency, which joins each deme to its immediate horizontal and vertical neighbors (the four cardinal directions), the same neighborhood a raster analysis with four directions would use. ```{r graph} coords <- as.matrix(expand.grid(x = 0:2, y = 0:1)) # 6 demes on a 3 x 2 grid rownames(coords) <- paste0("deme", seq_len(nrow(coords))) coords g <- deme_graph(coords, neighbours = "lattice") g ``` **How to read the output:** printing the graph reports the vertex count and the number of undirected edges. The object itself is a list with the parts you will pass downstream: ```{r graph-parts} str(g) g$edge_pairs ``` - `vertex_coordinates` is the `n x 2` coordinate matrix, one row per deme, in the order you supplied. - `edge_pairs` is an `m x 2` integer matrix of undirected edges. Each row holds the 1-based row indices of the two demes an edge joins. Each pair appears once, always with the smaller index first (`a < b`), and the rows are sorted. This canonical form is what makes the directional covariates later in this vignette unambiguous. - `n_vertices` is the deme count, a convenience copy of `nrow(vertex_coordinates)`. The returned object carries class `c("landgraph", "terradish_graph")`, so it drops directly into `dragonflow::dragon()` and is interchangeable with a `terradish::conductance_surface()` result. A quick plot makes the adjacency concrete. Edges are drawn first so the deme markers sit on top. ```{r graph-plot, fig.alt = "Six demes on a grid joined by rook-adjacency edges."} op <- par(no.readonly = TRUE) par(mar = c(4, 4, 1, 1)) plot(g$vertex_coordinates, type = "n", asp = 1, xlab = "x", ylab = "y") ep <- g$edge_pairs segments(g$vertex_coordinates[ep[, 1], 1], g$vertex_coordinates[ep[, 1], 2], g$vertex_coordinates[ep[, 2], 1], g$vertex_coordinates[ep[, 2], 2], col = "grey60") points(g$vertex_coordinates, pch = 21, bg = "steelblue", cex = 3) text(g$vertex_coordinates, labels = seq_len(nrow(coords)), col = "white") par(op) ``` **How to read this plot:** each blue circle is a deme, numbered by its row in `coords`. A grey line is an edge in `edge_pairs`. Rook adjacency gives the seven edges you see: the within-row horizontal links and the between-row vertical links, with no diagonals. ## Choosing the neighbourhood The `neighbours` argument selects how edges are drawn. The three options trade off assumptions about your sampling layout. | Value | What it does | When to use it | |-------|--------------|----------------| | `"lattice"` | Rook adjacency on a regular grid; add diagonals with `queen = TRUE` | Demes lie on an integer grid with uniform spacing | | `"knn"` | Joins each deme to its `k` nearest neighbors, then symmetrizes | Irregularly placed demes; `k` controls connectivity | | `"delaunay"` | The Delaunay triangulation (needs the `deldir` package) | Irregular demes, when you want a planar, parameter-free graph | For `"knn"`, set `k` to the number of neighbors each deme should reach; the result is symmetrized, so if deme A lists B among its neighbors the edge A-B is kept even when B does not list A. For `"lattice"`, leave `queen = FALSE` for the four cardinal neighbors, or set `queen = TRUE` to add the four diagonal neighbors (eight in total). `"delaunay"` takes no tuning but requires `deldir`; if the package is absent, `deme_graph()` stops with an instructive message rather than failing silently. ```{r graph-knn} g_knn <- deme_graph(coords, neighbours = "knn", k = 2) nrow(g_knn$edge_pairs) # edge count under 2-nearest-neighbour adjacency ``` # Genetic covariance from biallelic (SNP) data With a graph in hand, the next input is a genetic covariance matrix among the same demes. `cov_from_biallelic()` builds one from counts of the derived allele at biallelic (two-state) markers such as SNPs. We simulate counts for the six demes across eight SNPs, sampling 40 haploid chromosomes (20 diploid individuals) per deme. ```{r snp-data} set.seed(42) n_demes <- nrow(coords) n_snp <- 8 freqs <- runif(n_snp, 0.1, 0.9) # a true frequency per SNP Y <- vapply(freqs, function(p) rbinom(n_demes, size = 40, prob = p), numeric(n_demes)) rownames(Y) <- rownames(coords) colnames(Y) <- paste0("snp", seq_len(n_snp)) Y ``` `Y` is the **derived-allele count matrix**: demes in rows, loci in columns. Each cell is the number of copies of the derived (counted) allele observed in that deme at that locus. For diploid individuals genotyped 0/1/2, this is the standard allele-dosage matrix. The second input, `N`, is the **haploid sample size**: how many chromosomes were scored in each cell. You can supply it flexibly, and the choice you make should reflect how your sampling actually varied. | `N` form | Meaning | |----------|---------| | `NULL` (default) | Use `ploidy` for every cell (assumes complete genotyping) | | a single number | The same sample size in every cell | | length `nrow(Y)` | One size per deme, constant across loci | | length `ncol(Y)` | One size per locus, constant across demes | | a matrix matching `Y` | A separate size for every cell (handles missing data) | Here every cell was scored at 40 chromosomes, so a single number is enough. ```{r snp-cov} S <- cov_from_biallelic(Y, N = 40) round(S, 3) ``` **How to read the output:** `S` is a symmetric deme-by-deme covariance matrix on the scale of normalized allele frequencies. Internally the function standardizes each locus to the pooled allele frequency across demes, so each SNP contributes on a comparable scale, then averages the cross-products over loci (the genomic relationship matrix of Yang et al. 2010). Read the entries as relatedness in allele-frequency space: - A **diagonal** entry is a deme's variance, how far its standardized allele frequencies sit from the pooled mean. Larger means more distinctive. - An **off-diagonal** entry is the covariance between two demes. Positive means they deviate from the pooled frequencies in the same direction (genetically similar); negative means they deviate oppositely (differentiated). The matrix is positive semi-definite up to numerical tolerance, and row and column names are carried through from `Y`. This is exactly the response `S` that `terradish::wishart_covariance()` expects. A locus that is **monomorphic** (fixed at frequency 0 or 1 across all demes) carries no information and cannot be standardized. By default such loci are dropped with a warning; set `monomorphic = "error"` to stop instead. The `tol` argument sets how close to 0 or 1 a pooled frequency must be to count as fixed. # Pairwise F_ST For a more classical summary of differentiation, `fst_from_biallelic()` returns pairwise F_ST, the proportion of total genetic variation that is due to differences between demes rather than within them. It uses the ratio-of-averages estimator of Bhatia et al. (2013), which combines information across loci before taking the ratio. This function needs `N` as a full matrix. ```{r fst} Nmat <- matrix(40, n_demes, n_snp) fst <- fst_from_biallelic(Y, Nmat) round(fst, 4) ``` **How to read the output:** `fst` is a symmetric matrix with a zero diagonal (a deme has no differentiation from itself). Each off-diagonal entry is the pairwise F_ST between two demes: 0 means the pair is genetically indistinguishable, and larger positive values mean stronger differentiation. F_ST is already a proportion, so no back-transformation is needed. One caution on scale: the estimator is not constrained to `[0, 1]`. For very similar demes, sampling noise can push an estimate slightly below zero. Read a small negative value as "no detectable differentiation," not as an error. # Covariance from multivariate or microsatellite data Not all data are biallelic. `cov_from_genetic_data()` builds a covariance from any numeric genetic encoding: microsatellite allele calls, multiallelic markers, SNP dosages, or principal-component scores. It works at two levels. With no groups, each row is an individual and you get an individual-level covariance. With `groups`, rows are pooled to population centroids and you get a population-level covariance, the construction behind Dyer-style population graphs (Dyer and Nason 2004). The method is **Gower double-centering** (Gower 1966): it forms squared Euclidean distances among the units in feature space, then centers them into a covariance. We start from a small numeric feature matrix for three populations of two individuals each. ```{r multi-features} x <- matrix(c(0, 1, 1, 1, 2, 0, 2, 1, 0, 2, 1, 2), ncol = 2, byrow = TRUE) groups <- rep(c("pop1", "pop2", "pop3"), each = 2) Sg <- cov_from_genetic_data(x, groups = groups) round(Sg, 3) ``` **How to read the output:** `Sg` is a population-by-population covariance. Off the diagonal it behaves like the SNP covariance above: positive means two populations sit on the same side of the overall mean in feature space. The diagonal, however, is special here. Because these populations have replication (two individuals each), `cov_from_genetic_data()` defaults to `diagonal = "within"`, which replaces each diagonal entry with that population's **within-population genetic variance**, the spread of its members in feature space. This matches the covariance used before partial-correlation filtering in population-graph workflows. The function attaches the intermediate quantities as attributes, so you can inspect or reuse them. The most useful are listed below. ```{r multi-attrs} attr(Sg, "level") # "population" once groups have replication attr(Sg, "diagonal") # the diagonal rule actually applied attr(Sg, "within_variance") # the within-population variances on the diagonal attr(Sg, "unit_size") # number of individuals per population ``` | Attribute | What it holds | |-----------|---------------| | `level` | `"individual"` or `"population"` | | `diagonal` | The diagonal rule used: `"within"` or `"gower"` | | `within_variance` | Within-population variance per group (the `"within"` diagonal) | | `centroids` | Population centroids in feature space | | `centroid_distance2` | Squared distances among centroids | | `unit_size` | Individuals per group | | `retained_features` | Features kept after constant ones were dropped | The `diagonal` argument controls this behavior directly. `"auto"` (the default) uses `"within"` when groups have replication and `"gower"` otherwise. Force `"gower"` to keep the plain double-centered diagonal, or `"within"` to require the within-population variance. Two more arguments shape the feature space before centering: `center` and `scale` (both `TRUE` by default) standardize each feature, and `normalize = "features"` divides the result by the number of retained features so its scale does not grow with marker count. ## Microsatellite allele calls Microsatellite data usually arrive as allele calls, two columns per locus for a diploid. Set `input = "allele_calls"` and pass `loci` to tell the function which columns belong to the same locus. The calls are converted to per-allele dosage columns before centering. ```{r multi-msat} alleles <- data.frame( loc1_a = c(100, 100, 102, 102, 104, 104), loc1_b = c(100, 102, 102, 104, 104, 100), loc2_a = c(200, 202, 200, 202, 204, 204), loc2_b = c(202, 202, 204, 204, 204, 200) ) Sm <- cov_from_genetic_data( alleles, groups = groups, input = "allele_calls", loci = c("loc1", "loc1", "loc2", "loc2") ) round(Sm, 3) ``` **How to read the output:** the result is the same kind of population covariance as before. The `loci` vector has one entry per column of `alleles` and names the locus each allele copy belongs to; here the first two columns are `loc1` and the last two are `loc2`. A missing call is imputed to the most common (modal) allele observed at that locus, and the function reports how many calls it filled in. A modeling note for downstream Wishart fits. The effective degrees of freedom `nu` is the number of independent pieces of information in the covariance. For SNPs that is roughly the retained SNP count. For microsatellites it is safest to use the number of **loci**, because the allele frequencies within one locus are correlated (they sum to a constant) and so do not each count as independent. Report the value you use and check that conclusions hold across the plausible range. # From covariance to distance Some models want a distance matrix rather than a covariance. `dist_from_cov()` converts one to the other using the identity that relates a covariance to its implied squared Euclidean distances: $$D_{ij} = C_{ii} + C_{jj} - 2 C_{ij}.$$ ```{r dist} D <- dist_from_cov(S) round(D, 3) ``` **How to read the output:** `D` is a symmetric squared-distance matrix with a zero diagonal and non-negative off-diagonal entries. A larger value means two demes are farther apart in genetic space, the natural response for an isolation-by-distance or resistance model such as `terradish::mlpe()` or `terradish::generalized_wishart()`. Because `D` is built from `S`, the two describe the same structure: where the covariance is high, the distance is low. For biallelic data you can go straight from counts to distance with `dist_from_biallelic()`, a convenience wrapper for `dist_from_cov(cov_from_biallelic(Y, N))`. ```{r dist-wrap} D2 <- dist_from_biallelic(Y, N = 40) all.equal(D, D2) # same as the two-step route above ``` # Directional edge covariates The last piece is what makes directed gene-flow models possible. A **directional edge covariate** assigns a value to each edge that flips sign when you traverse the edge the other way, so it can describe an asymmetry such as flow downhill or with a prevailing wind. The covariate is **antisymmetric**: the value from deme `a` to deme `b` is the negative of the value from `b` to `a`. `landgraph` builds two kinds. ## Gradient of a scalar potential `edge_gradient()` takes a single value per deme (a **scalar potential** such as elevation) and returns the drop across each directed edge, `x_a - x_b`. Movement from high to low potential is the "downhill" direction. We use the sum of the coordinates as a stand-in elevation. ```{r grad} elevation <- coords[, "x"] + coords[, "y"] # one value per deme elevation eg <- edge_gradient(elevation, g) str(eg) ``` **How to read the output:** `edge_gradient()` returns a list with two parts. `edges` is an integer matrix of **directed** edges; every undirected edge from the graph appears twice, once in each direction (columns `a` and `b` are the start and end deme). `d` is the matching vector of potential drops, `x_a - x_b`, one per directed edge. Because each edge appears in both directions, the second half of `d` is the exact negative of the first half. That sign flip is the antisymmetry, and it is what a directed model uses to tell "uphill" from "downhill." A positive coefficient on this covariate in `terradish` directional models means movement speeds up as potential drops. ## Projection of a vector flow field `edge_gradient()` can only describe forces that point "downhill" from some potential. A real flow such as wind or current can also rotate, circling without any high or low point to descend from. `edge_flow()` captures that. It takes a **vector field** (an x-component and a y-component at each deme) and projects the average field along each edge onto the edge's direction. The covariate for undirected edge `(a, b)` is $$c_{ab} = \tfrac{1}{2}(f_a + f_b) \cdot (xy_b - xy_a),$$ the mean field on the edge dotted with the step from `a` to `b`. We build a counter-clockwise rotational field centered on the demes, the kind of pattern a scalar potential cannot represent. ```{r flow} cen <- colMeans(g$vertex_coordinates) rotational <- function(xy) cbind(-(xy[, 2] - cen[2]), xy[, 1] - cen[1]) circ <- edge_flow(rotational, g) round(circ, 3) ``` **How to read the output:** `circ` has one entry per undirected edge in `g$edge_pairs`, in the same row order. The sign tells you whether the field pushes along the edge from `a` to `b` (positive) or from `b` to `a` (negative), and the magnitude is how strongly. The covariate is antisymmetric by construction: a downstream model applies `circ` to the `a -> b` direction and its negative to `b -> a`. Pass it straight to `dragonflow::dragon()` as `circulation = circ`. The `field` argument also accepts a two-column matrix of components in vertex order, or a two-layer `terra::SpatRaster` sampled at the deme coordinates; the function form used here is convenient when you can express the field analytically. The plot below draws the field as an arrow at each deme, over the graph. The arrows circle the center, which is exactly the rotational signal `edge_flow()` extracts and `edge_gradient()` would miss. ```{r flow-plot, fig.alt = "A counter-clockwise vector field drawn as arrows at each deme over the graph edges."} op <- par(no.readonly = TRUE) par(mar = c(4, 4, 1, 1)) vc <- g$vertex_coordinates fld <- rotational(vc) plot(vc, type = "n", asp = 1, xlab = "x", ylab = "y", xlim = range(vc[, 1]) + c(-0.4, 0.4), ylim = range(vc[, 2]) + c(-0.4, 0.4)) segments(vc[ep[, 1], 1], vc[ep[, 1], 2], vc[ep[, 2], 1], vc[ep[, 2], 2], col = "grey80") arrows(vc[, 1], vc[, 2], vc[, 1] + 0.3 * fld[, 1], vc[, 2] + 0.3 * fld[, 2], length = 0.08, col = "firebrick") points(vc, pch = 21, bg = "steelblue", cex = 2.5) par(op) ``` **How to read this plot:** each red arrow is the flow vector at a deme; together they trace a counter-clockwise circulation. `edge_flow()` reads this field along each grey edge and returns the signed strength in `circ`. Where an arrow points along an edge, that edge gets a large covariate; where the flow crosses an edge sideways, the covariate is near zero. # Quick reference The complete pipeline, from coordinates and molecular data to the inputs a downstream model consumes: ```{r quick-ref, eval = FALSE} library(landgraph) # 1. Build the spatial graph from deme coordinates g <- deme_graph(coords, neighbours = "lattice") # or "knn" / "delaunay" # 2. Genetic covariance and distance from biallelic (SNP) counts S <- cov_from_biallelic(Y, N = 40) # deme covariance D <- dist_from_cov(S) # squared genetic distance D <- dist_from_biallelic(Y, N = 40) # the two steps in one call fst <- fst_from_biallelic(Y, Nmat) # pairwise F_ST (N as a matrix) # 2b. Covariance from microsatellite / multivariate data Sg <- cov_from_genetic_data(x, groups = groups) # population covariance Sg <- cov_from_genetic_data(alleles, groups = groups, input = "allele_calls", loci = c("loc1", "loc1", "loc2", "loc2")) # 3. Directional edge covariates for directed models eg <- edge_gradient(elevation, g) # downhill drop; eg$d per directed edge circ <- edge_flow(rotational, g) # circulation per undirected edge # 4. Hand off downstream # terradish::wishart_covariance(S = S, ...) # symmetric resistance # dragonflow::dragon(..., circulation = circ) # asymmetric gene flow ``` # Summary of key functions | Function | Purpose | |----------|---------| | `deme_graph()` | Build a deme/landscape graph (vertices + undirected edges) from coordinates | | `cov_from_biallelic()` | Genetic covariance from biallelic (SNP) allele counts | | `fst_from_biallelic()` | Pairwise F_ST from biallelic allele counts | | `cov_from_genetic_data()` | Covariance from multivariate or microsatellite data | | `dist_from_cov()` | Convert a covariance matrix to a squared-distance matrix | | `dist_from_biallelic()` | Covariance-to-distance shortcut for biallelic counts | | `edge_gradient()` | Directional covariate from the gradient of a scalar potential | | `edge_flow()` | Directional covariate from the projection of a vector flow field | # See also - The function help pages, for the full argument lists and return values: `?deme_graph`, `?cov_from_biallelic`, `?cov_from_genetic_data`, `?edge_gradient`, and `?edge_flow`. - `terradish`, for symmetric landscape resistance models that consume the covariance and distance matrices built here. - `dragonflow`, for asymmetric (directional) gene-flow models that consume the graph and the directional edge covariates from `edge_gradient()` and `edge_flow()`. # References Bhatia G, Patterson N, Sankararaman S, Price AL. 2013. Estimating and interpreting F_ST: the impact of rare variants. Genome Research 23(9):1514-1521. doi:10.1101/gr.154831.113 Dyer RJ, Nason JD. 2004. Population graphs: the graph theoretic shape of genetic structure. Molecular Ecology 13(7):1713-1727. doi:10.1111/j.1365-294X.2004.02177.x Gower JC. 1966. Some distance properties of latent root and vector methods used in multivariate analysis. Biometrika 53(3-4):325-338. doi:10.1093/biomet/53.3-4.325 Yang J, Benyamin B, McEvoy BP, et al. 2010. Common SNPs explain a large proportion of the heritability for human height. Nature Genetics 42(7):565-569. doi:10.1038/ng.608