7  Genetic Distance & Relatedness

This chapter covers inter-individual and inter-population distance estimation, pairwise relatedness, and spatial autocorrelation analysis.

library(gstudio)
library(dplyr)
data(arapat)

7.1 Genetic Distance

The genetic_distance() function computes distance matrices from genetic data. Some modes operate among individuals, others among populations (strata).

7.1.1 Available Distance Metrics

Mode	Level	Description
"AMOVA"	Individual	Inter-individual AMOVA distance
"Bray"	Individual	Proportion of shared alleles
"Euclidean"	Population	Euclidean frequency distance
"cGD"	Population	Conditional genetic distance
"pGD"	Population	Partitioned (directed) conditional genetic distance; asymmetric
"Nei"	Population	Nei’s genetic distance (1978)
"Dps"	Population	Shared allele distance (1 - Ps)
"Jaccard"	Population	Jaccard set dissimilarity

7.1.2 Individual-Level Distances

AMOVA distance computes inter-individual genetic distances:

# The first 20 individuals, for a manageable example
D_amova <- arapat |>
  slice(1:20) |>
  genetic_distance(mode = "AMOVA")

dim(D_amova)

[1] 20 20

D_amova[1:5, 1:5]

          [,1]     [,2]     [,3]     [,4]     [,5]
[1,]  0.000000 7.272727 9.090909 3.636364 10.90909
[2,]  7.272727 0.000000 1.000000 2.000000  5.00000
[3,]  9.090909 1.000000 0.000000 3.000000  2.00000
[4,]  3.636364 2.000000 3.000000 0.000000  5.00000
[5,] 10.909091 5.000000 2.000000 5.000000  0.00000

7.1.3 Population-Level Distances

Population-level distances require a stratum column:

D_nei <- arapat |>
  genetic_distance(stratum = "Species", mode = "Nei")
D_nei

               Cape  Mainland Peninsula
Cape      0.0000000 0.7267072 0.9233605
Mainland  0.7267072 0.0000000 0.6812596
Peninsula 0.9233605 0.6812596 0.0000000

D_cgd <- arapat |>
  genetic_distance(stratum = "Species", mode = "cGD")

Warning in popgraph(x = mv, groups = factor(as.character(x[[stratum]]))): 1
variables are collinear and being dropped from the discriminant rotation.

D_cgd

              Cape  Mainland Peninsula
Cape      0.000000  9.913083  7.136108
Mainland  9.913083  0.000000 16.490315
Peninsula 7.136108 16.490315  0.000000

Partitioned conditional genetic distance (pGD) splits each Population Graph edge’s cGD into two directed shares using the genetic-gravity neighbourhood weights, so the matrix is asymmetric: [i, j] is the distance from i to j, shorter in the direction of gene flow. The bandwidth multiplier is set with gamma (default 0.5); see vignette("genetic_gravity").

D_pgd <- arapat |>
  genetic_distance(stratum = "Species", mode = "pGD")

Warning in popgraph(x = mv, groups = factor(as.character(x[[stratum]]))): 1
variables are collinear and being dropped from the discriminant rotation.

D_pgd

              Cape Mainland Peninsula
Cape      0.000000 1.936519  3.214012
Mainland  7.976564 0.000000 11.190575
Peninsula 3.922096 4.155494  0.000000

D_euc <- arapat |>
  genetic_distance(stratum = "Species", mode = "Euclidean")

Multilocus estimates of Euclidean distance are assumed to be additive.

D_euc

              Cape Mainland Peninsula
Cape      0.000000 2.418731  2.494205
Mainland  2.418731 0.000000  2.140258
Peninsula 2.494205 2.140258  0.000000

7.2 Genetic Relatedness

The genetic_relatedness() function computes pairwise relatedness among individuals.

7.2.1 Available Relatedness Estimators

Mode	Description
"Nason"	Nason’s kinship coefficient (Fij)
"LynchRitland"	Lynch & Ritland estimator
"Ritland"	Ritland estimator

7.2.2 Basic Usage

The Nason estimator divides by the polymorphic index (genetic_diversity(mode = "Pe")), so monomorphic loci (where Pe = 0) will produce NaN. In the arapat dataset, LTRS and ATPS can be monomorphic in small subsets. To avoid this, restrict the data to polymorphic loci:

# Use loci that are polymorphic in this subset
R <- arapat |>
  slice(1:10) |>
  select(Species, WNT, EN, EF, ZMP, AML, MP20) |>
  genetic_relatedness(mode = "Nason")

Some of your loci are missing, Fij will treat these as loci with all alleles with likelihood equal to the population allele frequency.
Some of your loci are missing, Fij will treat these as loci with all alleles with likelihood equal to the population allele frequency.
Some of your loci are missing, Fij will treat these as loci with all alleles with likelihood equal to the population allele frequency.

dim(R)

[1] 10 10

round(R[1:5, 1:5], 4)

        [,1]    [,2]    [,3]    [,4]    [,5]
[1,]  1.0000 -0.5978 -0.8941 -0.4274 -1.0200
[2,] -0.5978  1.0000  0.6948  0.4308  0.2256
[3,] -0.8941  0.6948  1.0000  0.1345  1.1145
[4,] -0.4274  0.4308  0.1345  1.0000 -0.1038
[5,] -1.0200  0.2256  1.1145 -0.1038  1.0000

7.2.3 Single-Locus Relatedness

Relatedness at a single locus uses the same function, with loci =:

r_nason <- arapat |>
  slice(1:10) |>
  genetic_relatedness(loci = "EN", mode = "Nason")

round(r_nason[1:5, 1:5], 4)

        [,1]    [,2]    [,3]    [,4]    [,5]
[1,]  1.0000 -0.9953 -0.9953  0.0271  0.0271
[2,] -0.9953  1.0000  0.9220 -0.1153 -0.1153
[3,] -0.9953  0.9220  1.0000 -0.1153 -0.1153
[4,]  0.0271 -0.1153 -0.1153  1.0000  0.2329
[5,]  0.0271 -0.1153 -0.1153  0.2329  1.0000

7.3 Spatial Autocorrelation

The genetic_autocorrelation() function performs the spatial autocorrelation analysis of Smouse & Peakall (1999). It requires two square distance matrices — one physical, one genetic — and a vector of distance bin boundaries:

# Physical distance matrix
P <- arapat |>
  strata_distance(longitude = "Longitude", latitude = "Latitude")

# Genetic distance matrix
G <- arapat |>
  genetic_distance(mode = "AMOVA")

# Define distance bins
bins <- seq(0, max(P), length.out = 6)

# Run autocorrelation
result <- genetic_autocorrelation(P, G, bins = bins, perms = 99)
result

The output includes the spatial autocorrelation coefficient (R) for each distance class, the number of pairwise comparisons (N), and a p-value from the permutation test (P).

7.4 Strata Distances

The strata_coordinates() function extracts per-stratum centroid coordinates, and strata_distance() computes physical (geographic) distances between them:

coords <- arapat |>
  strata_coordinates(stratum = "Species", longitude = "Longitude", latitude = "Latitude")
coords

    Stratum Longitude Latitude
3 Peninsula -111.9226 26.18956
1      Cape -110.5805 24.25405
2  Mainland -109.6759 26.97435

coords |>
  strata_distance()

          Peninsula     Cape Mainland
Peninsula    0.0000 254.0525 239.8471
Cape       254.0525   0.0000 315.7868
Mainland   239.8471 315.7868   0.0000

7.5 Without Pipes

D_amova <- genetic_distance(arapat[1:20, ], mode = "AMOVA")
D_nei   <- genetic_distance(arapat, stratum = "Species", mode = "Nei")
D_pgd   <- genetic_distance(arapat, stratum = "Species", mode = "pGD", gamma = 0.5)
R       <- genetic_relatedness(arapat[1:10, c("Species", "WNT", "EN", "EF", "ZMP", "AML", "MP20")],
                               mode = "Nason")
coords  <- strata_coordinates(arapat, stratum = "Species")
strata_distance(coords)