Getting Started with landgraph
Bill Peterman
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.
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"), 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.
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")
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.
"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.
g_knn <- deme_graph(coords, neighbours = "knn", k = 2)
nrow(g_knn$edge_pairs) # edge count under 2-nearest-neighbour adjacency
#> [1] 7
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.
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.
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
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.
Nmat <- matrix(40, n_demes, n_snp)
fst <- fst_from_biallelic(Y, Nmat)
round(fst, 4)
#> [,1] [,2] [,3] [,4] [,5] [,6]
#> [1,] 0.0000 0.0013 0.0026 -0.0129 0.0291 -0.0155
#> [2,] 0.0013 0.0000 0.0150 0.0033 0.0051 -0.0039
#> [3,] 0.0026 0.0150 0.0000 0.0126 0.0161 0.0087
#> [4,] -0.0129 0.0033 0.0126 0.0000 0.0310 -0.0047
#> [5,] 0.0291 0.0051 0.0161 0.0310 0.0000 0.0259
#> [6,] -0.0155 -0.0039 0.0087 -0.0047 0.0259 0.0000
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.
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)
#> pop1 pop2 pop3
#> pop1 0.625 -0.429 0.067
#> pop2 -0.429 0.882 -1.605
#> pop3 0.067 -1.605 0.625
#> attr(,"centroids")
#> feature1 feature2
#> pop1 -0.559017 -0.2214037
#> pop2 1.118034 -0.8856149
#> pop3 -0.559017 1.1070186
#> attr(,"unit_features")
#> feature1 feature2
#> pop1 -0.559017 -0.2214037
#> pop2 1.118034 -0.8856149
#> pop3 -0.559017 1.1070186
#> attr(,"within_variance")
#> pop1 pop2 pop3
#> 0.6250000 0.8823529 0.6250000
#> attr(,"centroid_distance2")
#> pop1 pop2 pop3
#> pop1 0.000000 3.253676 1.764706
#> pop2 3.253676 0.000000 6.783088
#> pop3 1.764706 6.783088 0.000000
#> attr(,"unit_distance2")
#> pop1 pop2 pop3
#> pop1 0.000000 3.253676 1.764706
#> pop2 3.253676 0.000000 6.783088
#> pop3 1.764706 6.783088 0.000000
#> attr(,"unit_size")
#> pop1 pop2 pop3
#> 2 2 2
#> attr(,"feature_center")
#> feature1 feature2
#> 1.000000 1.166667
#> attr(,"feature_scale")
#> feature1 feature2
#> 0.8944272 0.7527727
#> attr(,"retained_features")
#> [1] "feature1" "feature2"
#> attr(,"input")
#> [1] "features"
#> attr(,"level")
#> [1] "population"
#> attr(,"diagonal")
#> [1] "within"
#> attr(,"normalize")
#> [1] "none"
#> attr(,"normalizer")
#> [1] 1
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.
attr(Sg, "level") # "population" once groups have replication
#> [1] "population"
attr(Sg, "diagonal") # the diagonal rule actually applied
#> [1] "within"
attr(Sg, "within_variance") # the within-population variances on the diagonal
#> pop1 pop2 pop3
#> 0.6250000 0.8823529 0.6250000
attr(Sg, "unit_size") # number of individuals per population
#> pop1 pop2 pop3
#> 2 2 2
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.
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 3.917 -1.328 -2.689
#> pop2 -1.328 3.917 -0.512
#> pop3 -2.689 -0.512 4.049
#> attr(,"centroids")
#> loc1:100 loc1:102 loc1:104 loc2:200 loc2:202 loc2:204
#> pop1 1.0206207 -0.2041241 -0.8164966 0 1.0206207 -1.1070186
#> pop2 -0.8164966 1.0206207 -0.2041241 0 -0.2041241 0.2214037
#> pop3 -0.2041241 -0.8164966 1.0206207 0 -0.8164966 0.8856149
#> attr(,"unit_features")
#> loc1:100 loc1:102 loc1:104 loc2:200 loc2:202 loc2:204
#> pop1 1.0206207 -0.2041241 -0.8164966 0 1.0206207 -1.1070186
#> pop2 -0.8164966 1.0206207 -0.2041241 0 -0.2041241 0.2214037
#> pop3 -0.2041241 -0.8164966 1.0206207 0 -0.8164966 0.8856149
#> attr(,"within_variance")
#> pop1 pop2 pop3
#> 3.916667 3.916667 4.049020
#> attr(,"centroid_distance2")
#> pop1 pop2 pop3
#> pop1 0.000000 8.514706 12.595588
#> pop2 8.514706 0.000000 6.066176
#> pop3 12.595588 6.066176 0.000000
#> attr(,"unit_distance2")
#> pop1 pop2 pop3
#> pop1 0.000000 8.514706 12.595588
#> pop2 8.514706 0.000000 6.066176
#> pop3 12.595588 6.066176 0.000000
#> attr(,"unit_size")
#> pop1 pop2 pop3
#> 2 2 2
#> attr(,"feature_center")
#> loc1:100 loc1:102 loc1:104 loc2:200 loc2:202 loc2:204
#> 0.6666667 0.6666667 0.6666667 0.5000000 0.6666667 0.8333333
#> attr(,"feature_scale")
#> loc1:100 loc1:102 loc1:104 loc2:200 loc2:202 loc2:204
#> 0.8164966 0.8164966 0.8164966 0.5477226 0.8164966 0.7527727
#> attr(,"retained_features")
#> [1] "loc1:100" "loc1:102" "loc1:104" "loc2:200" "loc2:202" "loc2:204"
#> attr(,"input")
#> [1] "allele_calls"
#> attr(,"imputed_allele_calls")
#> named integer(0)
#> attr(,"imputed_modal_alleles")
#> named character(0)
#> attr(,"level")
#> [1] "population"
#> attr(,"diagonal")
#> [1] "within"
#> attr(,"normalize")
#> [1] "none"
#> attr(,"normalizer")
#> [1] 1
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}.\]
D <- dist_from_cov(S)
round(D, 3)
#> deme1 deme2 deme3 deme4 deme5 deme6
#> [1,] 0.000 1.996 2.122 0.984 3.615 0.806
#> [2,] 1.996 0.000 2.977 2.483 2.282 1.722
#> [3,] 2.122 2.977 0.000 3.291 3.144 2.532
#> [4,] 0.984 2.483 3.291 0.000 4.175 1.521
#> [5,] 3.615 2.282 3.144 4.175 0.000 3.671
#> [6,] 0.806 1.722 2.532 1.521 3.671 0.000
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)).
D2 <- dist_from_biallelic(Y, N = 40)
all.equal(D, D2) # same as the two-step route above
#> [1] TRUE
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.
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 ...
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.
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
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.
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)
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:
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
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