05 · Phylo-Geographic Analysis

Isolation by Distance (Mantel Test), Scatter Pie Map, and Mantel Correlogram

Performs isolation-by-distance analysis using Mantel tests on patristic distances from a phylogenetic tree and geographic (Haversine) distances. Produces a scatter pie map showing clade composition per sampling location and a Mantel correlogram.

Overview

Item	Details
Input	Aligned FASTA (S segment required, M optional) · Newick tree · coordinates CSV · metadata CSV
Key packages	ape, vegan, geosphere, scatterpie, sf, rnaturalearth, ggspatial
Statistics	Mantel test (Pearson) · patristic distances · Haversine geographic distances
Output	Isolation-by-distance plot · scatter pie map · Mantel correlogram · TIFF + SVG
Download	template.Rmd

Tip

When to use this template: You have a phylogenetic tree and sampling coordinates for a pathogen, and want to test whether genetic distance correlates with geographic distance (isolation by distance) and visualize the geographic distribution of clades.

Note

Clade colors: Default palette uses Okabe-Ito colorblind-safe colors. Update CLADE_COLORS in the USER CONFIGURATION block to match your clade assignments. The scatter pie radius (radius column in pie_data) is in degrees — adjust mutate(radius = 0.4) to scale pies appropriately for your map extent.

Back to Gallery Open Template File

## ── USER CONFIGURATION ──────────────────────────────────────────────────────
#
# Input files:
S_FASTA    <- "data/S_aligned.fasta"   # Aligned S-segment FASTA
M_FASTA    <- "data/M_aligned.fasta"   # Aligned M-segment FASTA  (optional;
                                       #   set to NULL to skip it entirely)
TREEFILE   <- "data/S_tree.nwk"        # Newick tree from IQ-TREE or similar
COORDS_FILE<- "data/coordinates.csv"   # Columns: SampleID, Longitude, Latitude
METADATA_FILE <- "data/metadata.csv"   # Columns: SampleID, Clade (+ optional extras)

# Column names in COORDS_FILE and METADATA_FILE that link to FASTA tip labels.
SAMPLE_COL <- "SampleID"

# TIP_ID_REGEX: How to recover the SampleID from a tree tip label.
#   IQ-TREE / RAxML tips routinely carry decoration, e.g.
#   "Sample_12|Clade_A|2019" or "Sample_12_Chile_2019". This regex is applied
#   to BOTH the tip labels and the metadata IDs, so the two sides always match
#   whether or not the decoration is present.
#   Set to NULL to use tip labels verbatim.
TIP_ID_REGEX <- "Sample_[0-9]+"

# Columns in METADATA_FILE used for clade coloring.
CLADE_COL  <- "Clade"

# CLADE_COLORS: Named color vector (Okabe-Ito palette).
#   Names must match the unique values in METADATA_FILE[[CLADE_COL]].
#   Clades found in the data but missing here are NOT dropped — they are given
#   fallback colours and reported, so the pies always sum to the true sample
#   composition.
CLADE_COLORS <- c(
  "Clade_A" = "#E69F00",   # Orange
  "Clade_B" = "#56B4E9",   # Sky blue
  "Clade_C" = "#009E73",   # Bluish green
  "Clade_D" = "#0072B2"    # Blue
)

# CLADE_FALLBACK_COLORS: Used, in order, for clades absent from CLADE_COLORS.
CLADE_FALLBACK_COLORS <- c("#D55E00", "#CC79A7", "#F0E442", "#999999", "#661100")

# PATHOGEN_NAME: Used in plot titles and axis labels.
PATHOGEN_NAME <- "Pathogen X"

# Mantel test: number of permutations (increase to 99999 for publication).
MANTEL_NPERM <- 999

# Map extent (degrees): lon/lat range around your sampling area.
MAP_LON <- c(-72, -63)
MAP_LAT <- c(-30, -20)

# OUTPUT_DIR: Where to save plots.
OUTPUT_DIR <- "Plots"

## ────────────────────────────────────────────────────────────────────────────

Setup

library(ape)
library(vegan)
library(geosphere)
library(Biostrings)
library(tidyverse)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(ggspatial)
library(scatterpie)
library(svglite)
library(patchwork)
library(ragg)

dir.create(file.path(OUTPUT_DIR, "Maps"),    recursive = TRUE, showWarnings = FALSE)
dir.create(file.path(OUTPUT_DIR, "Mantel"),  recursive = TRUE, showWarnings = FALSE)

Load Data

# FASTA sequences
seqs_S <- readDNAStringSet(S_FASTA)

# Phylogenetic tree
tree_S <- read.tree(TREEFILE)

# Metadata & coordinates
coords   <- read_csv(COORDS_FILE, show_col_types = FALSE)
metadata <- read_csv(METADATA_FILE, show_col_types = FALSE)

# Combine
sample_data <- coords %>%
  left_join(metadata, by = SAMPLE_COL) %>%
  filter(!is.na(Longitude), !is.na(Latitude))

cat("Sequences loaded:", length(seqs_S), "\n")

Sequences loaded: 30 

cat("Samples with coordinates:", nrow(sample_data), "\n")

Samples with coordinates: 30 

# Column names are taken from the data, so the documented "+ optional extras"
# in METADATA_FILE works instead of erroring on a hardcoded col.names length.
sample_data %>%
  head(10) %>%
  knitr::kable(
    caption   = "Sampling Coordinates and Metadata (First 10 Samples)",
    col.names = gsub("_", " ", names(sample_data))
  )

Sampling Coordinates and Metadata (First 10 Samples)

SampleID	Longitude	Latitude	Clade	Year	Host
Sample_1	-69.8	-21.5	Clade_A	2017	Rodent_B
Sample_2	-69.8	-21.5	Clade_A	2021	Rodent_B
Sample_3	-69.8	-21.5	Clade_A	2015	Rodent_B
Sample_4	-69.8	-21.5	Clade_A	2016	Rodent_A
Sample_5	-68.2	-22.8	Clade_A	2020	Rodent_A
Sample_6	-68.2	-22.8	Clade_A	2017	Rodent_A
Sample_7	-68.2	-22.8	Clade_A	2015	Rodent_A
Sample_8	-68.2	-22.8	Clade_B	2023	Rodent_A
Sample_9	-70.1	-24.0	Clade_B	2018	Rodent_B
Sample_10	-70.1	-24.0	Clade_B	2021	Rodent_B

Alignment Quality

# IUPAC degeneracy codes. MAFFT and several other aligners emit lowercase
# sequence, so the comparison is done on an upper-cased copy — matching only
# uppercase silently reports zero ambiguities for every such sequence.
IUPAC_AMBIG <- c("N", "R", "Y", "S", "W", "K", "M", "B", "D", "H", "V")

ambiguity_table <- function(seqs, segment_label) {
  chars <- strsplit(toupper(as.character(seqs)), "", fixed = TRUE)
  data.frame(
    Segment     = segment_label,
    Sequence    = names(seqs),
    Length_bp   = as.integer(width(seqs)),
    N_Bases     = vapply(chars, function(s) sum(s == "N"), integer(1)),
    Ambig_Bases = vapply(chars, function(s) sum(s %in% IUPAC_AMBIG), integer(1)),
    Gap_Bases   = vapply(chars, function(s) sum(s == "-"), integer(1)),
    row.names   = NULL,
    stringsAsFactors = FALSE
  )
}

qc_seq <- ambiguity_table(seqs_S, "S")

# The M segment is optional: include it in the same QC table when configured.
seqs_M <- NULL
if (!is.null(M_FASTA) && nzchar(M_FASTA) && file.exists(M_FASTA)) {
  seqs_M <- readDNAStringSet(M_FASTA)
  qc_seq <- rbind(qc_seq, ambiguity_table(seqs_M, "M"))
  cat("M-segment sequences loaded:", length(seqs_M), "\n")
} else {
  cat("M segment not analysed (M_FASTA is unset or missing).\n")
}

M-segment sequences loaded: 30 

qc_seq %>%
  mutate(Ambig_Pct = 100 * Ambig_Bases / Length_bp) %>%
  arrange(desc(Ambig_Pct)) %>%
  head(10) %>%
  knitr::kable(
    caption   = "Sequence Alignment Quality (10 Most Ambiguous Sequences)",
    col.names = c("Segment", "Sequence Header", "Length (bp)", "N Bases",
                  "Ambiguous Bases (all IUPAC codes)", "Gaps", "Ambiguous %"),
    digits    = 2
  )

Sequence Alignment Quality (10 Most Ambiguous Sequences)

Segment	Sequence Header	Length (bp)	N Bases	Ambiguous Bases (all IUPAC codes)	Gaps	Ambiguous %
S	Sample_14	1200	11	32	0	2.67
M	Sample_7	3000	20	78	0	2.60
M	Sample_28	3000	25	78	0	2.60
M	Sample_21	3000	16	76	0	2.53
S	Sample_21	1200	7	28	0	2.33
M	Sample_14	3000	14	70	0	2.33
S	Sample_28	1200	3	26	0	2.17
S	Sample_7	1200	4	24	0	2.00
M	Sample_3	3000	1	19	0	0.63
S	Sample_12	1200	3	7	0	0.58

qc_seq %>%
  group_by(Segment) %>%
  summarise(
    Sequences      = n(),
    `Alignment bp` = paste(unique(Length_bp), collapse = ", "),
    `Total ambiguous` = sum(Ambig_Bases),
    `Mean ambiguous %` = mean(100 * Ambig_Bases / Length_bp),
    .groups = "drop"
  ) %>%
  knitr::kable(caption = "Alignment Summary by Segment", digits = 2)

Alignment Summary by Segment

Segment	Sequences	Alignment bp	Total ambiguous	Mean ambiguous %
M	30	3000	443	0.49
S	30	1200	158	0.44

Isolation by Distance (Mantel Test)

# ── 0. Normalise IDs on BOTH sides ────────────────────────────────────────────
# The patristic matrix is named by raw tip labels; the geographic matrix by
# metadata SampleIDs. Applying the extraction to only one side means that the
# moment tips carry any decoration ("Sample_12|Clade_A|2019") the two name sets
# stop intersecting, the matrices collapse to 0x0, and mantel() receives empty
# input. Applying it to both sides is idempotent when there is no decoration.
extract_sample_id <- function(x) {
  if (is.null(TIP_ID_REGEX)) return(as.character(x))
  out <- str_extract(as.character(x), TIP_ID_REGEX)
  # Anything the regex does not match keeps its original label rather than
  # becoming NA, so a partial mismatch is visible instead of silent.
  ifelse(is.na(out), as.character(x), out)
}

# ── 1. Patristic distance matrix from tree ─────────────────────────────────────
patristic_mat <- cophenetic.phylo(tree_S)
rownames(patristic_mat) <- colnames(patristic_mat) <-
  extract_sample_id(rownames(patristic_mat))

# ── 2. Geographic distance matrix (great-circle, km) ──────────────────────────
sample_data <- sample_data %>%
  mutate(.sample_id = extract_sample_id(.data[[SAMPLE_COL]]))

# Order the coordinates to follow the tree tips, keeping only samples that have
# both a tip and a coordinate.
tip_sample_ids <- extract_sample_id(tree_S$tip.label)

coords_ordered <- sample_data %>%
  filter(.sample_id %in% tip_sample_ids) %>%
  distinct(.sample_id, .keep_all = TRUE) %>%
  slice(match(intersect(tip_sample_ids, .sample_id), .sample_id))

geo_mat <- distm(
  cbind(coords_ordered$Longitude, coords_ordered$Latitude),
  fun = distHaversine
) / 1000   # → km

rownames(geo_mat) <- colnames(geo_mat) <- coords_ordered$.sample_id

# ── 3. Mantel test ─────────────────────────────────────────────────────────────
shared  <- intersect(rownames(patristic_mat), rownames(geo_mat))

if (length(shared) < 3) {
  message("Tip labels: ", paste(head(tree_S$tip.label, 3), collapse = ", "), " ...")
  message("Metadata IDs: ", paste(head(sample_data[[SAMPLE_COL]], 3), collapse = ", "), " ...")
  message("TIP_ID_REGEX = ", if (is.null(TIP_ID_REGEX)) "NULL" else TIP_ID_REGEX)
}
stopifnot(
  "Fewer than 3 samples are shared between the tree and the coordinates - check TIP_ID_REGEX" =
    length(shared) >= 3
)

n_tree_only <- length(setdiff(rownames(patristic_mat), shared))
n_geo_only  <- length(setdiff(rownames(geo_mat), shared))
cat(sprintf("Mantel test on %d samples (%d tree-only, %d coordinate-only excluded)\n",
            length(shared), n_tree_only, n_geo_only))

Mantel test on 30 samples (0 tree-only, 0 coordinate-only excluded)

pat_sub <- patristic_mat[shared, shared]
geo_sub <- geo_mat[shared, shared]

mantel_result <- mantel(pat_sub, geo_sub,
                         method  = "pearson",
                         permutations = MANTEL_NPERM)

mantel_df <- tibble(
  `Mantel Statistic (r)` = unname(mantel_result$statistic),
  `p-value (Significance)` = unname(mantel_result$signif),
  `Permutations`          = MANTEL_NPERM,
  `Correlation Method`    = "Pearson"
)

mantel_df %>%
  knitr::kable(
    caption = paste("Mantel Test Result (Isolation by Distance) —", PATHOGEN_NAME),
    digits  = 4
  )

Mantel Test Result (Isolation by Distance) — Pathogen X

Mantel Statistic (r)	p-value (Significance)	Permutations	Correlation Method
0.8118	0.001	999	Pearson

# ── 4. Scatter plot ────────────────────────────────────────────────────────────
ibd_df <- data.frame(
  geographic  = as.vector(geo_sub[lower.tri(geo_sub)]),
  phylogenetic = as.vector(pat_sub[lower.tri(pat_sub)])
)

p_ibd <- ggplot(ibd_df, aes(x = geographic, y = phylogenetic)) +
  geom_point(alpha = 0.5, size = 1.5, color = "steelblue") +
  geom_smooth(method = "lm", color = "steelblue", se = TRUE) +
  labs(
    title = paste("Isolation by Distance —", PATHOGEN_NAME),
    subtitle = sprintf("Mantel r = %.3f,  p = %.3f (n_perm = %d)",
                       mantel_result$statistic, mantel_result$signif, MANTEL_NPERM),
    x = "Geographic distance (km)",
    y = "Patristic distance (branch length)"
  ) +
  theme_bw(base_size = 10, base_family = "sans") +
  theme(
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

ggsave(file.path(OUTPUT_DIR, "Mantel", "isolation_by_distance.tiff"),
       plot = p_ibd, width = 3, height = 3, dpi = 300,
       device = "tiff", compression = "lzw")
print(p_ibd)

Scatter Pie Map

# ── Summarize clade composition per location ───────────────────────────────────
# (when multiple samples share a location, show their clade proportions as a pie)
pie_data <- sample_data %>%
  group_by(Longitude, Latitude, !!sym(CLADE_COL)) %>%
  summarise(n = n(), .groups = "drop") %>%
  tidyr::pivot_wider(names_from = all_of(CLADE_COL), values_from = n, values_fill = 0) %>%
  mutate(radius = 0.4)   # pie radius in degrees; adjust to your map scale

# ── Palette completion ────────────────────────────────────────────────────────
# Plotting only intersect(names(CLADE_COLORS), names(pie_data)) drops any clade
# the config forgot to list. The pies then still look correct while no longer
# representing the true sample composition. Instead: keep every clade in the
# data and give the unlisted ones fallback colours, loudly.
clades_present <- sort(unique(as.character(sample_data[[CLADE_COL]])))
clade_palette  <- setNames(unname(CLADE_COLORS[clades_present]), clades_present)

missing_clades <- clades_present[is.na(clade_palette)]
if (length(missing_clades) > 0) {
  clade_palette[missing_clades] <-
    rep(CLADE_FALLBACK_COLORS, length.out = length(missing_clades))
  message("Clades present in the data but absent from CLADE_COLORS (assigned ",
          "fallback colours): ", paste(missing_clades, collapse = ", "))
}

unused_clades <- setdiff(names(CLADE_COLORS), clades_present)
if (length(unused_clades) > 0) {
  message("CLADE_COLORS entries not present in the data (ignored): ",
          paste(unused_clades, collapse = ", "))
}

pie_cols <- clades_present
stopifnot(all(pie_cols %in% names(pie_data)))

# Every sample must be represented by exactly one pie slice.
stopifnot(sum(as.matrix(pie_data[, pie_cols])) == nrow(sample_data))

# Country polygons
world <- ne_countries(scale = "medium", returnclass = "sf")

p_spie <- ggplot() +
  geom_sf(data = world, fill = "grey90", color = "white", linewidth = 0.2) +
  geom_scatterpie(
    data    = pie_data,
    aes(x = Longitude, y = Latitude, r = radius),
    cols    = pie_cols,
    alpha   = 0.85, color = "black", linewidth = 0.2
  ) +
  scale_fill_manual(values = clade_palette, name = "Clade") +
  coord_sf(xlim = MAP_LON, ylim = MAP_LAT, expand = FALSE) +
  annotation_scale(location = "br", width_hint = 0.25) +
  labs(title = paste(PATHOGEN_NAME, "— Geographic Distribution"),
       x = NULL, y = NULL) +
  theme_bw(base_size = 9, base_family = "sans") +
  theme(
    legend.position   = "top",
    legend.direction  = "horizontal",
    legend.text       = element_text(size = 6),
    legend.key.size   = unit(3, "pt"),
    panel.grid.major  = element_blank(),
    panel.grid.minor  = element_blank(),
    panel.border      = element_rect(color = "black", linewidth = 0.35),
    axis.text         = element_text(size = 8),
    plot.margin       = margin(4, 4, 4, 4, "pt")
  )

ggsave(file.path(OUTPUT_DIR, "Maps", "scatter_pie.tiff"),
       plot = p_spie, width = 3, height = 2, dpi = 600,
       device = "tiff", compression = "lzw")
ggsave(file.path(OUTPUT_DIR, "Maps", "scatter_pie.svg"),
       plot = p_spie, width = 3, height = 2, device = "svg")
print(p_spie)

Mantel Correlogram

# Distance class boundaries (km). Quartiles of the observed pairwise distances
# give four classes with equal numbers of pairs — a reasonable default for any
# study area. Replace `probs` (or supply explicit break points) to match yours.
dist_breaks <- quantile(as.vector(geo_sub[lower.tri(geo_sub)]),
                         probs = seq(0, 1, by = 0.25))

# Mantel correlogram (vegan). break.pts is what makes the classes above take
# effect; without it vegan picks its own Sturges-rule classes and dist_breaks
# is computed for nothing.
mcor <- mantel.correlog(
  D.eco    = pat_sub,
  D.geo    = geo_sub,
  break.pts = as.numeric(dist_breaks),
  nperm    = MANTEL_NPERM,
  cutoff   = FALSE
)

mcor_df <- as.data.frame(mcor$mantel.res, check.names = FALSE)

# vegan reports raw and, when mult != "none", corrected p-values.
# Match the package's own plotting behavior by preferring corrected p-values.
p_col <- if ("Pr(corrected)" %in% names(mcor_df)) "Pr(corrected)" else "Pr(Mantel)"

mcor_df <- mcor_df %>%
  mutate(
    pvalue = .data[[p_col]],
    sig = ifelse(!is.na(pvalue) & pvalue <= 0.05, "Significant", "NS")
  )

p_mcor <- ggplot(mcor_df, aes(x = `class.index`, y = Mantel.cor, fill = sig)) +
  geom_bar(stat = "identity", width = diff(range(mcor_df$`class.index`)) / nrow(mcor_df) * 0.8) +
  geom_hline(yintercept = 0, color = "black", linewidth = 0.4) +
  scale_fill_manual(values = c("Significant" = "steelblue", "NS" = "grey70"),
                    name = NULL) +
  labs(title = paste("Mantel Correlogram —", PATHOGEN_NAME),
       x = "Geographic distance class (km)",
       y = "Mantel r") +
  theme_bw(base_size = 10, base_family = "sans") +
  theme(
    legend.position = c(0.8, 0.75),
    legend.background = element_rect(color = "black"),
    legend.key = element_rect(fill = "transparent"),
    panel.grid.major = element_blank(),
    panel.grid.minor = element_blank()
  )

ggsave(file.path(OUTPUT_DIR, "Mantel", "mantel_correlogram.tiff"),
       plot = p_mcor, width = 3, height = 2, dpi = 300,
       device = "tiff", compression = "lzw")
print(p_mcor)