1 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 terradish, which estimates symmetric landscape resistance. 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.

library(landgraph)

2 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.

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
#>       x y
#> deme1 0 0
#> deme2 1 0
#> deme3 2 0
#> deme4 0 1
#> deme5 1 1
#> deme6 2 1

g <- deme_graph(coords, neighbours = "lattice")
g
#> landgraph: 6 vertices, 7 undirected edges

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:

str(g)
#> List of 4
#>  $ vertex_coordinates: num [1:6, 1:2] 0 1 2 0 1 2 0 0 0 1 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:6] "deme1" "deme2" "deme3" "deme4" ...
#>   .. ..$ : chr [1:2] "x" "y"
#>  $ edge_pairs        : int [1:7, 1:2] 1 1 2 2 3 4 5 2 4 3 ...
#>  $ coords            : num [1:6, 1:2] 0 1 2 0 1 2 0 0 0 1 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : chr [1:6] "deme1" "deme2" "deme3" "deme4" ...
#>   .. ..$ : chr [1:2] "x" "y"
#>  $ n_vertices        : int 6
#>  - attr(*, "class")= chr [1:2] "landgraph" "terradish_graph"
g$edge_pairs
#>      [,1] [,2]
#> [1,]    1    2
#> [2,]    1    4
#> [3,]    2    3
#> [4,]    2    5
#> [5,]    3    6
#> [6,]    4    5
#> [7,]    5    6

• 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"). The class lets edge_gradient() and edge_flow() accept the object. It lacks the covariates, Laplacian, and Cholesky factor that terradish::terradish() needs; build fitting graphs with terradish::conductance_surface().

A quick plot makes the adjacency concrete. Edges are drawn first so the deme markers sit on top.

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 = "gray60")
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 gray 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.

g_knn <- deme_graph(coords, neighbours = "knn", k = 2)
nrow(g_knn$edge_pairs)   # edge count under 2-nearest-neighbor adjacency
#> [1] 7

3 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.

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
#>       snp1 snp2 snp3 snp4 snp5 snp6 snp7 snp8
#> deme1   32   34   17   31   25   16   28    7
#> deme2   32   30   10   27   23   20   22    7
#> deme3   34   29   20   31   32   22   24    8
#> deme4   32   37   18   28   22   21   27   10
#> deme5   30   34    9   29   32   26   22    4
#> deme6   35   34   13   28   27   15   27   10

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.

Here every cell was scored at 40 chromosomes, so a single number is enough.

S <- cov_from_biallelic(Y, N = 40)
round(S, 3)
#>        deme1  deme2  deme3  deme4  deme5  deme6
#> deme1  0.551 -0.286 -0.132  0.303 -0.643  0.208
#> deme2 -0.286  0.873 -0.398 -0.285  0.184 -0.088
#> deme3 -0.132 -0.398  1.308 -0.472 -0.029 -0.276
#> deme4  0.303 -0.285 -0.472  1.039 -0.679  0.095
#> deme5 -0.643  0.184 -0.029 -0.679  1.778 -0.610
#> deme6  0.208 -0.088 -0.276  0.095 -0.610  0.672
#> attr(,"diagonal")
#> [1] "normalized_dosage"
#> attr(,"centered")
#> [1] "sites"

How to read the output: S is a symmetric deme-by-deme covariance matrix of normalized allele counts. 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. The two attributes printed below the matrix record how it was built. The function standardizes counts, not frequencies, so the entries grow with the haploid sample size. Every deme here was scored at 40 chromosomes, so S / 40 gives the standardized allele-frequency covariance that terradish fits. Do not divide by a common N when sample sizes differ; the demes then carry unequal sampling variance that a common nugget does not absorb.

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.

4 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. N takes the same forms as in cov_from_biallelic(), so a single number works here.

fst <- fst_from_biallelic(Y, N = 40)
round(fst, 4)
#>         deme1   deme2  deme3   deme4  deme5   deme6
#> deme1  0.0000  0.0013 0.0026 -0.0129 0.0291 -0.0155
#> deme2  0.0013  0.0000 0.0150  0.0033 0.0051 -0.0039
#> deme3  0.0026  0.0150 0.0000  0.0126 0.0161  0.0087
#> deme4 -0.0129  0.0033 0.0126  0.0000 0.0310 -0.0047
#> deme5  0.0291  0.0051 0.0161  0.0310 0.0000  0.0259
#> deme6 -0.0155 -0.0039 0.0087 -0.0047 0.0259  0.0000

How to read the output: fst is a symmetric matrix named by the demes, 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.

5 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.

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)   # Sg[, ] drops the attributes for a clean printout
#>        pop1   pop2   pop3
#> pop1  0.362 -0.429  0.067
#> pop2 -0.429  2.034 -1.605
#> pop3  0.067 -1.605  1.538

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. Because cov_from_genetic_data() defaults to diagonal = "gower", the diagonal holds each population’s variance on the same scale as the covariances, rows sum to zero, and the matrix is positive semidefinite, as a Wishart likelihood requires. Population-graph workflows instead use diagonal = "within", whose diagonal is on another scale and does not suit Wishart models.

The function attaches the intermediate quantities as attributes, so you can inspect or reuse them. The most useful are listed below.

attr(Sg, "level")            # "population" once groups have replication
#> [1] "population"
attr(Sg, "diagonal")         # the diagonal rule actually applied
#> [1] "gower"
attr(Sg, "within_variance")  # fills the diagonal only with diagonal = "within"
#>      pop1      pop2      pop3 
#> 0.6250000 0.8823529 0.6250000
attr(Sg, "unit_size")        # number of individuals per population
#> pop1 pop2 pop3 
#>    2    2    2

The diagonal argument controls this behavior directly. "gower" (the default) keeps the double-centered diagonal; "within" substitutes within-population variance; "auto" restores the default used before 0.0.3. Two more arguments shape the feature space before centering. center and scale (both TRUE by default) standardize each feature. Centering leaves the covariance unchanged, while scaling gives rare alleles the same variance as common ones and so more relative weight. normalize = "features" divides the result by the number of retained allele features so its scale does not grow with marker count; that divisor is not nu.

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)
#>        pop1   pop2   pop3
#> pop1  4.017 -1.328 -2.689
#> pop2 -1.328  1.841 -0.512
#> pop3 -2.689 -0.512  3.201

6 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}.\]

D <- dist_from_cov(S)
round(D, 3)
#>       deme1 deme2 deme3 deme4 deme5 deme6
#> deme1 0.000 1.996 2.122 0.984 3.615 0.806
#> deme2 1.996 0.000 2.977 2.483 2.282 1.722
#> deme3 2.122 2.977 0.000 3.291 3.144 2.532
#> deme4 0.984 2.483 3.291 0.000 4.175 1.521
#> deme5 3.615 2.282 3.144 4.175 0.000 3.671
#> deme6 0.806 1.722 2.532 1.521 3.671 0.000
#> attr(,"diagonal")
#> [1] "normalized_dosage"
#> attr(,"centered")
#> [1] "sites"

How to read the output: D is a symmetric squared-distance matrix with a zero diagonal and non-negative off-diagonal entries, named by the demes. It keeps the diagonal and centered attributes of S. 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. D inherits the count scale of S, so with 40 chromosomes per deme the terradish response dist_from_cov(S / 40) equals D / 40.

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)).

D2 <- dist_from_biallelic(Y, N = 40)
all.equal(D, D2)                     # same as the two-step route above
#> [1] TRUE
all.equal(dist_from_cov(S / 40), D / 40)   # frequency scale for terradish
#> [1] TRUE

7 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.

elevation <- coords[, "x"] + coords[, "y"]   # one value per deme
elevation
#> deme1 deme2 deme3 deme4 deme5 deme6 
#>     0     1     2     1     2     3

eg <- edge_gradient(elevation, g)
str(eg)
#> List of 2
#>  $ edges: int [1:14, 1:2] 1 1 2 2 3 4 5 2 4 3 ...
#>   ..- attr(*, "dimnames")=List of 2
#>   .. ..$ : NULL
#>   .. ..$ : chr [1:2] "a" "b"
#>  $ d    : num [1:14] -1 -1 -1 -1 -1 -1 -1 1 1 1 ...

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)
#> deme1 deme1 deme2 deme2 deme3 deme4 deme5 
#>   0.5  -1.0   0.5   0.0   1.0  -0.5  -0.5

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 = "gray80")
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)

8 Quick reference

The complete pipeline, from coordinates and molecular data to the inputs a downstream model consumes:

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, N = 40)               # pairwise F_ST (any N form)

# 2b. Covariance from microsatellite / multivariate data
Sg <- cov_from_genetic_data(x, groups = groups)    # population covariance, Gower diagonal
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 (equal N: divide by N for the frequency scale)
#    terradish fits the response S / 40 (or dist_from_cov(S / 40))
#    with nu supplied by you; see terradish::terradish_rescale_nu()

9 Summary of key functions

10 See also

• The function help pages, for the full argument lists and return values: ?deme_graph, ?cov_from_biallelic, ?fst_from_biallelic, ?cov_from_genetic_data, ?dist_from_cov, ?dist_from_biallelic, ?edge_gradient, and ?edge_flow.
• terradish, for symmetric landscape resistance models that consume the covariance and distance matrices built here; start with vignette("getting-started", package = "terradish").

11 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