Beau Larkin
Last updated: 04 October, 2026
- Description
- Sequence data processing functions
- Confidence interval helper
- Select spatial eigenvectors
- Alpha diversity calculations
- Confidence intervals
- Multivariate analysis
- Permanova on soil data
- Model distribution probabilities
- Filter spe to a guild
- Calculate pairwise distances among sites and present summary statistics
- Site mapping functions
Functions that accompany the repository, sourced from this file to save space elsewhere
etl <- function(spe, env, taxa, traits = NULL, varname, gene, cluster_type = "otu",
colname_prefix, folder) {
varname <- enquo(varname)
data <- spe %>% left_join(taxa, by = join_by(`#OTU ID`))
meta <- if (gene == "ITS") {
data %>%
mutate(!!varname := paste0(cluster_type, "_", row_number())) %>%
select(!starts_with(gene)) %>%
rename(otu_ID = `#OTU ID`) %>%
select(!!varname, everything()) %>%
separate_wider_delim(taxonomy, delim = ";",
names = c("kingdom", "phylum", "class", "order", "family", "genus", "species"),
too_few = "align_start") %>%
mutate(across(kingdom:species, ~ str_sub(.x, 4))) %>%
left_join(traits, by = join_by(phylum, class, order, family, genus)) %>%
select(-kingdom, -Confidence)
} else {
data %>%
mutate(!!varname := paste0(cluster_type, "_", row_number())) %>%
select(!starts_with(gene)) %>%
rename(otu_ID = `#OTU ID`) %>%
select(!!varname, everything()) %>%
separate(taxonomy,
into = c("class", "order", "family", "genus", "taxon", "accession"),
sep = ";", fill = "right") %>%
select(-Confidence)
}
spe_samps <- data %>%
mutate(!!varname := paste0(cluster_type, "_", row_number())) %>%
select(!!varname, starts_with(gene)) %>%
column_to_rownames(var = as_name(varname)) %>%
t() %>% as.data.frame() %>% rownames_to_column("rowname") %>%
mutate(rowname = str_remove(rowname, colname_prefix)) %>%
separate_wider_delim("rowname", delim = "_", names = c("field_key", "sample")) %>%
mutate(across(c(field_key, sample), as.numeric)) %>%
left_join(env %>% select(field_key, field_name), by = join_by(field_key)) %>%
select(field_name, sample, everything(), -field_key) %>%
arrange(field_name, sample)
spe_avg <- spe_samps %>%
group_by(field_name) %>%
summarize(across(starts_with(cluster_type), mean), .groups = "drop") %>%
arrange(field_name)
samples_fields <- spe_samps %>%
count(field_name, name = "n") %>%
left_join(sites, by = join_by(field_name)) %>%
select(field_name, region, n) %>%
arrange(n, field_name) %>%
kable(format = "pandoc", caption = paste("Number of samples per field", gene, sep = ",\n"))
write_csv(meta, root_path(folder, paste("spe", gene, "metadata.csv", sep = "_")))
write_csv(spe_samps, root_path(folder, paste("spe", gene, "samples.csv", sep = "_")))
write_csv(spe_avg, root_path(folder, paste("spe", gene, "avg.csv", sep = "_")))
list(samples_fields = samples_fields, spe_meta = meta, spe_samps = spe_samps, spe_avg = spe_avg)
}spe_accum <- function(data) {
df <- data.frame(
samples = specaccum(data[, -c(1, 2)], method = "exact")$site,
richness = specaccum(data[, -c(1, 2)], method = "exact")$richness,
sd = specaccum(data[, -c(1, 2)], method = "exact")$sd
)
df
}ci <- function(x) std.error(x) * qnorm(0.975)Fit null and full db-RDA models using dbMEM spatial variables and forward-select spatial eigenvectors associated with community composition.
mem_select <- function(d, mem, seed = 20260211, permutations = 1999) {
stopifnot(identical(labels(d), rownames(mem)))
mod_null <- dbrda(d ~ 1, data = mem)
mod_full <- dbrda(d ~ ., data = mem)
set.seed(seed)
ordistep(
mod_null,
scope = formula(mod_full),
direction = "forward",
permutations = permutations,
trace = FALSE
)
}Returns a dataframe of alpha diversity (richness, Shannon’s) for analysis and plotting. Handles the biofuel plot collapse internally
calc_div <- function(spe, site_dat, biofuel_plots = c("FLRSP1", "FLRSP2", "FLRSP3")) {
# Calculate sequencing depth and alpha diversity for each sampled plot
div_data <- spe %>%
rowwise() %>%
mutate(
depth = sum(c_across(starts_with("otu"))),
richness = sum(c_across(starts_with("otu")) > 0),
shannon = exp(diversity(c_across(starts_with("otu"))))
) %>%
select(field_name, depth, richness, shannon) %>%
ungroup()
# Retain ordinary sites unchanged
div_other <- div_data %>%
filter(!field_name %in% biofuel_plots) %>%
mutate(
depth_rich = depth,
depth_shan = depth
) %>%
select(field_name, depth_rich, depth_shan, richness, shannon)
# Collapse Fermi biofuel control plots to one independent replicate
biofuel <- div_data %>%
filter(field_name %in% biofuel_plots)
if (nrow(biofuel) > 0) {
median_richness <- median(biofuel$richness)
div_biofuel <- tibble(
field_name = "FLRSP1",
depth_rich = biofuel %>%
filter(richness == median_richness) %>%
summarize(depth = mean(depth)) %>%
pull(depth),
depth_shan = mean(biofuel$depth),
richness = median_richness,
shannon = mean(biofuel$shannon)
)
div_data <- bind_rows(div_other, div_biofuel)
} else {
div_data <- div_other
}
# Join site metadata and prepare transformed sequencing-depth covariates
div_data %>%
left_join(
site_dat %>% select(field_type, field_name),
by = join_by(field_name)
) %>%
mutate(
depth_rich_csq = sqrt(depth_rich) - mean(sqrt(depth_rich)),
depth_shan_csq = sqrt(depth_shan) - mean(sqrt(depth_shan))
) %>%
select(
field_name,
field_type,
depth_rich,
depth_rich_csq,
depth_shan,
depth_shan_csq,
richness,
shannon
)
}Calculate upper and lower confidence intervals with alpha=0.05
ci_u <- function(x) {(sd(x) / sqrt(length(x))) * qnorm(0.975)}
ci_l <- function(x) {(sd(x) / sqrt(length(x))) * qnorm(0.025)}NMDS ordination → dispersion check → global & pairwise PERMANOVA Args: d dist, env metadata, covar optional covariates (MEM), nperm permutations.
mva <- function(d, env, covar = NULL, nperm = 1999, seed = 20260211, plot_stress = TRUE) {
stopifnot(is.data.frame(env))
if (!("field_type" %in% names(env))) stop("`env` must contain column `field_type`.")
if (is.matrix(d)) {
if (!isTRUE(all.equal(d, t(d)))) stop("`d` matrix must be symmetric.")
diag(d) <- 0
}
# covariate checks
if (!is.null(covar)) {
if (!is.character(covar)) stop("`covar` must be NULL or a character vector of column names.")
covar <- unique(covar)
if (length(covar) < 1L || length(covar) > 3L) stop("`covar` must be NULL or 1–3 column name strings.")
missing_cov <- setdiff(covar, names(env))
if (length(missing_cov) > 0L) stop("`env` is missing covariate column(s): ", paste(missing_cov, collapse = ", "))
if (anyNA(env[, covar, drop = FALSE])) {
bad <- covar[colSums(is.na(env[, covar, drop = FALSE])) > 0]
stop("Covariate column(s) contain NA: ", paste(bad, collapse = ", "), "; handle before calling mva().")
}
}
# Distance labels and env alignment
if (inherits(d, "dist")) {
lab <- attr(d, "Labels")
} else if (is.matrix(d)) {
if (is.null(rownames(d))) stop("Distance matrix `d` must have row names.")
lab <- rownames(d)
d <- as.dist(d)
} else {
stop("`d` must be a 'dist' or a symmetric distance matrix.")
}
# Align feature names across data sources
env <- as.data.frame(env)
if ("field_name" %in% names(env)) rownames(env) <- env$field_name
if (!setequal(rownames(env), lab)) {
miss_env <- setdiff(lab, rownames(env))
miss_d <- setdiff(rownames(env), lab)
stop("Sample mismatch between `d` and `env`.\n",
"In d not in env: ", paste(miss_env, collapse = ", "),
"\nIn env not in d: ", paste(miss_d, collapse = ", "))
}
env <- env[lab, , drop = FALSE]
env$field_name <- rownames(env)
# Set grouping variables
g_chr <- as.character(env$field_type)
g_levels <- sort(unique(g_chr))
clust_vec <- factor(g_chr, levels = g_levels)
# Ordination (NMDS)
if (!is.null(seed)) set.seed(seed + 1L)
p <- metaMDS(
d,
k = 2,
trymax = 100,
autotransform = FALSE,
trace = FALSE
)
p_sco <- scores(
p,
display = "sites",
choices = 1:2
) %>%
as.data.frame() %>%
rownames_to_column(var = "field_name") %>%
left_join(env, by = join_by(field_name))
# Homogeneity of multivariate dispersion
disper <- betadisper(d, clust_vec, bias.adjust = TRUE)
if (!is.null(seed)) set.seed(seed + 2L)
mvdisper <- permutest(disper, pairwise = TRUE, permutations = nperm)
# Global PERMANOVA
if (!is.null(seed)) set.seed(seed + 3L)
perm_terms <- c(covar, "field_type")
perm_form <- reformulate(perm_terms, response = "d")
gl_permtest <- adonis2(
perm_form,
data = env,
permutations = nperm,
by = "terms"
)
# Pairwise PERMANOVA
groups <- combn(g_levels, m = 2) %>% t() %>% as.data.frame()
names(groups) <- c("V1", "V2")
contrasts <- data.frame(
group1 = groups$V1,
group2 = groups$V2,
R2 = NA_real_,
F_value = NA_real_,
df1 = NA_integer_,
df2 = NA_integer_,
p_value = NA_real_
)
d_mat <- as.matrix(d)
for (i in seq_len(nrow(contrasts))) {
g1 <- contrasts$group1[i]
g2 <- contrasts$group2[i]
keep <- clust_vec %in% c(g1, g2)
contrast_mat <- d_mat[keep, keep, drop = FALSE]
env_sub <- env[keep, , drop = FALSE]
env_sub$field_type <- droplevels(factor(env_sub$field_type))
if (!is.null(seed)) set.seed(seed + 100L + i)
perm_terms_pw <- c(covar, "field_type")
perm_form_pw <- reformulate(perm_terms_pw, response = "contrast_mat")
fit <- adonis2(
perm_form_pw,
data = env_sub,
permutations = nperm,
by = "terms"
)
rn <- rownames(fit)
term_row <- grep("(^field_type$)|field_type", rn)
if (length(term_row) != 1L) {
stop("Could not uniquely identify `field_type` row in pairwise adonis2 result.\nRows were: ",
paste(rn, collapse = ", "))
}
contrasts$R2[i] <- round(fit[term_row, "R2"], 3)
contrasts$F_value[i] <- round(fit[term_row, "F"], 3)
contrasts$df1[i] <- fit[term_row, "Df"]
contrasts$df2[i] <- fit[grep("^Residual", rn), "Df"]
contrasts$p_value[i] <- fit[term_row, "Pr(>F)"]
}
contrasts$p_value_adj <- round(p.adjust(contrasts$p_value, method = "fdr"), 4)
par(mfrow = c(1,1))
if (plot_stress) stressplot(p)
list(
ordination = p,
stress = p$stress,
ordination_scores = p_sco,
dispersion_test = mvdisper,
permanova = gl_permtest,
pairwise_contrasts = contrasts
)
}Simplified version of mva() for use with the soil properties data
soilperm <- function(d, env, covar = NULL, nperm = 1999, seed = 20251103) {
# Distance labels and env alignment
if (inherits(d, "dist")) {
lab <- attr(d, "Labels")
} else if (is.matrix(d)) {
if (is.null(rownames(d))) stop("Distance matrix `d` must have row names.")
lab <- rownames(d)
d <- as.dist(d) # coerce for betadisper/pcoa convenience
} else {
stop("`d` must be a 'dist' or a symmetric distance matrix.")
}
# Align feature names across data sources
env <- as.data.frame(env)
if ("field_name" %in% names(env)) rownames(env) <- env$field_name
if (!setequal(rownames(env), lab)) {
miss_env <- setdiff(lab, rownames(env))
miss_d <- setdiff(rownames(env), lab)
stop("Sample mismatch between `d` and `env`.\n",
"In d not in env: ", paste(miss_env, collapse = ", "),
"\nIn env not in d: ", paste(miss_d, collapse = ", "))
}
env <- env[lab, , drop = FALSE]
# Set grouping variables
g_chr <- as.character(env$field_type)
g_levels <- sort(unique(g_chr)) # deterministic order
clust_vec <- factor(g_chr, levels = g_levels)
# Homogeneity of multivariate dispersion
disper <- betadisper(d, clust_vec, bias.adjust = TRUE)
if (!is.null(seed)) set.seed(seed + 2L)
mvdisper <- permutest(disper, pairwise = TRUE, permutations = nperm)
# Global PERMANOVA
if (!is.null(seed)) set.seed(seed + 3L)
perm_terms <- c(covar, "field_type")
perm_form <- reformulate(perm_terms, response = "d")
gl_permtest <- adonis2(
perm_form,
data = env,
permutations = nperm,
by = "terms"
)
# Pairwise PERMANOVA
groups <- combn(g_levels, m = 2) %>% t() %>% as.data.frame()
names(groups) <- c("V1", "V2")
contrasts <- data.frame(
group1 = groups$V1,
group2 = groups$V2,
R2 = NA_real_,
F_value = NA_real_,
df1 = NA_integer_,
df2 = NA_integer_,
p_value = NA_real_
)
d_mat <- as.matrix(d)
for (i in seq_len(nrow(contrasts))) {
g1 <- contrasts$group1[i]
g2 <- contrasts$group2[i]
keep <- clust_vec %in% c(g1, g2)
contrast_mat <- d_mat[keep, keep, drop = FALSE]
env_sub <- env[keep, , drop = FALSE]
env_sub$field_type <- droplevels(factor(env_sub$field_type))
if (!is.null(seed)) set.seed(seed + 100L + i)
perm_terms_pw <- c(covar, "field_type")
perm_form_pw <- reformulate(perm_terms_pw, response = "contrast_mat")
fit <- adonis2(
perm_form_pw,
data = env_sub,
permutations = nperm,
by = "terms"
)
rn <- rownames(fit)
term_row <- grep("(^field_type$)|field_type", rn)
if (length(term_row) != 1L) {
stop("Could not uniquely identify `field_type` row in pairwise adonis2 result.\nRows were: ",
paste(rn, collapse = ", "))
}
contrasts$R2[i] <- round(fit[term_row, "R2"], 3)
contrasts$F_value[i] <- round(fit[term_row, "F"], 3)
contrasts$df1[i] <- fit[term_row, "Df"]
contrasts$df2[i] <- fit[grep("^Residual", rn), "Df"]
contrasts$p_value[i] <- fit[term_row, 5]
}
contrasts$p_value_adj <- p.adjust(contrasts$p_value, method = "fdr") %>% round(4)
# Results
list(
dispersion_test = mvdisper,
permanova = gl_permtest,
pairwise_contrasts = contrasts
)
}Probable distributions of response and residuals. Package performance prints javascript which doesn’t render on github documents.
distribution_prob <- function(df) {
print(
performance::check_distribution(df) %>%
as.data.frame() %>%
select(Distribution, p_Residuals) %>%
arrange(-p_Residuals) %>%
slice_head(n = 3) %>%
kable(format = "pandoc"))
print(
performance::check_distribution(df) %>%
as.data.frame() %>%
select(Distribution, p_Response) %>%
arrange(-p_Response) %>%
slice_head(n = 3) %>%
kable(format = "pandoc"))
}Create samp-spe matrix of sequence abundance in a guild
guildseq <- function(spe, meta, guild) {
guab <-
spe %>%
pivot_longer(starts_with("otu"), names_to = "otu_num", values_to = "abund") %>%
left_join(meta %>% select(otu_num, primary_lifestyle), by = join_by(otu_num)) %>%
filter(primary_lifestyle == guild) %>%
select(-primary_lifestyle) %>%
pivot_wider(names_from = otu_num, values_from = abund)
return(guab)
}Function reg_dist_stats() requires a haversine distance matrix, site
metadata, and is filtered by regions to produce the desired output.
reg_dist_stats <- function(dist_mat,
sites_df,
filt_rg) {
dm <- round(as.matrix(dist_mat) / 1000, 1)
if (!setequal(rownames(dm), as.character(sites_df$field_key))) {
stop("Distance matrix labels do not match sites_df$field_key")
}
idx <- which(upper.tri(dm), arr.ind = TRUE)
pairs <- tibble(
site1 = rownames(dm)[idx[, 1]],
site2 = colnames(dm)[idx[, 2]],
dist = dm[idx]
)
meta <- sites_df %>%
select(site = field_key, ft = field_type, rg = region) %>%
mutate(
site = as.character(site),
ft = as.character(ft)
) %>%
filter(rg == filt_rg)
pairs %>%
filter(site1 %in% meta$site, site2 %in% meta$site) %>%
left_join(
meta %>% select(site, ft),
by = c("site1" = "site")
) %>%
rename(ft1 = ft) %>%
left_join(
meta %>% select(site, ft),
by = c("site2" = "site")
) %>%
rename(ft2 = ft) %>%
mutate(
group_pair = paste(
pmin(ft1, ft2),
pmax(ft1, ft2),
sep = "-"
)
) %>%
group_by(group_pair) %>%
summarize(
min_dist = min(dist, na.rm = TRUE),
median_dist = median(dist, na.rm = TRUE),
max_dist = max(dist, na.rm = TRUE),
.groups = "drop"
)
}Functions facilitate the creation of mapping objects
bbox_buffer_km <- function(pts_sf, buffer_km = 20) {
albers <- 5070 # NAD83 / Conus Albers
pts_sf %>%
st_transform(albers) %>%
st_bbox() %>%
st_as_sfc(crs = albers) %>%
st_buffer(dist = buffer_km * 1000) %>%
st_transform(4326) %>%
st_bbox()
}Add field names manually: - nudge_e_m = meters to move EAST (positive = east, negative = west) - nudge_n_m = meters to move NORTH (positive = north, negative = south)
nudge_coords <- function(sites, lon_col = "long", lat_col = "lat",
east_col = "nudge_e_m", north_col = "nudge_n_m") {
stopifnot(all(c(lon_col, lat_col) %in% names(sites)))
# If nudge columns are missing, treat as zeros
if (!(east_col %in% names(sites))) sites[[east_col]] <- 0
if (!(north_col %in% names(sites))) sites[[north_col]] <- 0
# Vectorized meters → degrees (WGS84 approximation)
m_per_deg_lat <- 111320 # ~ meters per degree latitude
m_per_deg_lon <- m_per_deg_lat * cospi(sites[[lat_col]] / 180)
dx_deg <- sites[[east_col]] / m_per_deg_lon
dy_deg <- sites[[north_col]] / m_per_deg_lat
sites %>%
mutate(
long_plot = .data[[lon_col]] + ifelse(is.finite(dx_deg), dx_deg, 0),
lat_plot = .data[[lat_col]] + ifelse(is.finite(dy_deg), dy_deg, 0)
)
}get_osm_roads <- function(bb, density = 8) {
# validate density
types <- c("motorway","trunk","primary","secondary",
"tertiary","unclassified","residential","service")
if (!is.numeric(density) || length(density) != 1L || !is.finite(density)) {
stop("`density` must be a single finite number in 1:8")
}
density <- as.integer(max(1L, min(length(types), density)))
vals <- types[seq_len(density)]
# build layer
res <-
opq(bbox = bb) %>%
add_osm_feature(key = "highway", value = vals) %>%
osmdata_sf()
return(st_crop(res$osm_lines, st_as_sfc(bb)))
}make_zoom_map <- function(bb, panel_tag = NULL, pos = c(0,1), show_counties = FALSE, road_data = NULL) {
crop_states <- st_crop(cont, bb)
crop_counties <- if (show_counties) st_crop(st_transform(counties, 4326), bb) else NULL
roads <- road_data
pts <- sites_sf %>%
filter(long >= bb["xmin"], long <= bb["xmax"],
lat >= bb["ymin"], lat <= bb["ymax"]) %>%
st_drop_geometry()
# Labels only where yr_since is available (restored sites)
pts_lab <- sites_plot %>%
filter(!is.na(yr_since)) %>%
mutate(lbl = as.character(round(yr_since, 0)))
g <- ggplot() +
geom_sf(data = crop_states, fill = "ivory", color = "black", linewidth = 0.5) +
{ if (!is.null(crop_counties)) geom_sf(data = crop_counties, fill = NA, color = "gray85", linewidth = 0.3) } +
{ if (!is.null(roads)) geom_sf(data = roads, color = "grey70", linewidth = 0.3) } +
geom_point(
data = sites_plot,
aes(x = long_plot, y = lat_plot, fill = field_type),
shape = 21, size = sm_size, stroke = lw, color = "black"
) +
geom_text(
data = pts_lab,
aes(x = long_plot, y = lat_plot, label = lbl),
size = yrtx_size, family = "sans", fontface = 2, color = "black"
) +
scale_fill_manual(values = ft_pal) +
annotation_scale(location = "bl", width_hint = 0.35, height = grid::unit(0.15, "cm")) +
coord_sf(
xlim = c(bb["xmin"], bb["xmax"]),
ylim = c(bb["ymin"], bb["ymax"]),
expand = FALSE
) +
labs(tag = panel_tag) +
theme_void() +
theme(
panel.background = element_rect(fill = "aliceblue", color = "black", linewidth = 0.5),
legend.position = "none",
plot.tag = element_text(size = 14, face = 1, hjust = 0),
plot.tag.position = pos
)
g
}add_fig7_rug <- function(p, comp_df,
y0, h,
forb_fill, grass_fill, v_nudge=0) {
comp_df <- comp_df %>%
mutate(ymin = y0 + v_nudge,
ymid = (y0 + forb_comp * h) + v_nudge,
ymax = (y0 + h) + v_nudge)
p +
geom_ribbon(
data = comp_df,
aes(x = gf_axis, ymin = ymin, ymax = ymid),
inherit.aes = FALSE,
fill = forb_fill
) +
geom_ribbon(
data = comp_df,
aes(x = gf_axis, ymin = ymid, ymax = ymax),
inherit.aes = FALSE,
fill = grass_fill
) +
coord_cartesian(clip = "off")
}map_pt_to_mm <- function(pt) pt / ggplot2::.pt