Modeling site, platform, relatedness, and spatial structure
Selçuk Korkmaz
2026-11-05
Beyond the classic subject / batch / study / time relations, splitGraph models several further leakage axes, in two families:
- Cluster-style relations — collection site, tissue / anatomical region, sequencing platform, and assay — are categorical groupings. Each is auto-detected from a metadata column and handled by its own constraint mode, exactly like
subject or batch.
- Pairwise relations — genetic relatedness and spatial proximity — are continuous and defined between pairs. They are modelled as thresholded edges and grouped by transitive closure, a partition a single categorical column cannot express.
This vignette builds and groups by each, shows how the threshold drives the pairwise grouping, and points out the two ways a structure-aware grouping goes wrong — collapsing into one group, or dissolving into singletons.
Cluster-style relations: site, region, platform, assay
graph_from_metadata() auto-detects site_id, region_id, platform_id, and assay_id columns and builds the corresponding typed nodes and edges. Each then has its own constraint mode, and all four behave identically — only the column and the mode name change.
meta <- data.frame(
sample_id = paste0("S", 1:6),
subject_id = c("P1", "P1", "P2", "P2", "P3", "P3"),
site_id = c("NYC", "NYC", "BOS", "BOS", "NYC", "BOS"),
region_id = c("cortex", "cortex", "cortex",
"hippocampus", "hippocampus", "hippocampus"),
platform_id = c("illumina", "illumina", "nanopore",
"nanopore", "illumina", "nanopore"),
assay_id = c("rnaseq", "rnaseq", "rnaseq", "wgs", "wgs", "wgs"),
stringsAsFactors = FALSE
)
g <- graph_from_metadata(meta, graph_name = "structure-demo")
cluster_modes <- c("site", "region", "platform", "assay")
do.call(cbind, lapply(
stats::setNames(cluster_modes, cluster_modes),
function(m) grouping_vector(derive_split_constraints(g, mode = m))
))
#> site region platform assay
#> S1 "site:NYC" "region:cortex" "platform:illumina" "assay:rnaseq"
#> S2 "site:NYC" "region:cortex" "platform:illumina" "assay:rnaseq"
#> S3 "site:BOS" "region:cortex" "platform:nanopore" "assay:rnaseq"
#> S4 "site:BOS" "region:hippocampus" "platform:nanopore" "assay:wgs"
#> S5 "site:NYC" "region:hippocampus" "platform:illumina" "assay:wgs"
#> S6 "site:BOS" "region:hippocampus" "platform:nanopore" "assay:wgs"
Each column is a different, equally defensible partition of the same six samples. Choosing between them is a scientific question, not a technical one: which of these axes must a model generalise across?
Site structure also has a validation rule of its own. A subject whose samples were collected at more than one site is flagged, because grouping by site alone would then place one individual on both sides of a split:
report <- validate_graph(g)
report$issues[report$issues$code == "subject_cross_site_overlap",
c("severity", "message")]
#> severity message
#> 4 warning Subject `subject:P3` has samples collected at multiple sites.
# Which constraint modes actually sever that path:
risks <- summarize_leakage_risks(g, constraint = derive_split_constraints(g, "subject"))
as.data.frame(risks)[as.data.frame(risks)$category == "subject_cross_site_overlap",
c("category", "severed")]
#> category severed
#> 4 subject_cross_site_overlap TRUE
Whatever mode is primary, every detected cluster relation is also carried into the split_spec as a blocking annotation, so a downstream consumer can block on site, region, platform, or assay even when the split unit is something else — here, subject:
spec <- as_split_spec(derive_split_constraints(g, mode = "subject"), graph = g)
spec$block_vars
#> [1] "site_group" "region_group" "platform_group" "assay_group"
head(spec$sample_data[, c("sample_id", "group_id", "site_group",
"region_group", "platform_group", "assay_group")])
#> sample_id group_id site_group region_group platform_group assay_group
#> 1 S1 subject:P1 NYC cortex illumina rnaseq
#> 2 S2 subject:P1 NYC cortex illumina rnaseq
#> 3 S3 subject:P2 BOS cortex nanopore rnaseq
#> 4 S4 subject:P2 BOS hippocampus nanopore wgs
#> 5 S5 subject:P3 NYC hippocampus illumina wgs
#> 6 S6 subject:P3 BOS hippocampus nanopore wgs
Any of these relations can also participate in a composite derivation, where several dependency sources are combined and each connected component becomes one group. Site and platform partition this cohort the same way, so grouping on both at once still leaves two groups:
constraint <- derive_split_constraints(
g, mode = "composite", strategy = "strict",
via = c("Site", "Platform")
)
grouping_vector(constraint)
#> S1 S2 S3 S4 S5
#> "component_1" "component_1" "component_2" "component_2" "component_1"
#> S6
#> "component_2"
Watch for a composite that collapses
A strict composite can only ever merge groups, never split them, so every relation you add risks merging everything. Add subject to the same call and the entire cohort becomes one group:
collapsed <- derive_split_constraints(
g, mode = "composite", strategy = "strict",
via = c("Site", "Platform", "Subject")
)
collapsed$metadata$n_groups
#> [1] 1
grouping_vector(collapsed)
#> S1 S2 S3 S4 S5
#> "component_1" "component_1" "component_1" "component_1" "component_1"
#> S6
#> "component_1"
The validation report predicted this. subject_cross_site_overlap fired above because subject P3 has one sample at each site; in a composite closure that subject is a bridge, so the NYC and BOS groups fuse and there is nothing left to hold out. One group is not an error — it is an arithmetically correct answer to the question that was asked — and splitGraph does not warn about it, so read metadata$n_groups before you rely on a composite.
Two ways out, both already on this page: leave the bridging relation out of via (the Site + Platform call above), or keep the coarse axis as a blocking annotation rather than the split unit, which is what site_group in the spec is for. vignette("faq-design-notes") works through the trade-off.
Pairwise relation: genetic relatedness
Some leakage is pairwise and continuous rather than a clean grouping. Genetic relatedness is the canonical example: a kinship coefficient — typically from a tool such as KING or PLINK — links pairs of subjects. relatedness_edges_from_kinship() takes such a pair table, keeps pairs at or above a threshold, and emits subject_related_to edges; mode = "relatedness" then groups by transitive closure over those edges (so a chain of related individuals lands in one group).
# A kinship table over subject pairs (one sample per subject here for clarity).
# P1-P2 and P2-P3 clear the threshold and chain together; P5-P6 form a second
# related pair; P1-P4 is too weak to count.
kin <- data.frame(
id1 = c("P1", "P2", "P1", "P5"),
id2 = c("P2", "P3", "P4", "P6"),
kinship = c(0.25, 0.20, 0.02, 0.30),
stringsAsFactors = FALSE
)
rel_edges <- relatedness_edges_from_kinship(kin, threshold = 0.1)
meta_r <- data.frame(
sample_id = paste0("S", 1:6),
subject_id = paste0("P", 1:6),
stringsAsFactors = FALSE
)
samples <- create_nodes(meta_r, "Sample", "sample_id")
subjects <- create_nodes(meta_r, "Subject", "subject_id")
belongs <- create_edges(meta_r, "sample_id", "subject_id",
"Sample", "Subject", "sample_belongs_to_subject")
g_rel <- build_dependency_graph(list(samples, subjects), list(belongs, rel_edges))
rel_groups <- grouping_vector(derive_split_constraints(g_rel, mode = "relatedness"))
rel_groups
#> S1 S2 S3
#> "relatedness:component_1" "relatedness:component_1" "relatedness:component_1"
#> S4 S5 S6
#> "relatedness:component_2" "relatedness:component_3" "relatedness:component_3"
Real kinship tools rarely use those exact column names, and they do not all emit the long format. Both shapes are accepted. KING and GCTA write a long table with their own headers, which id1 / id2 / kinship rename:
king <- data.frame(
ID1 = c("P1", "P2", "P1", "P5"),
ID2 = c("P2", "P3", "P4", "P6"),
Kinship = c(0.25, 0.20, 0.02, 0.30),
stringsAsFactors = FALSE
)
king_edges <- relatedness_edges_from_kinship(
king, threshold = 0.1, id1 = "ID1", id2 = "ID2", kinship = "Kinship"
)
as.data.frame(king_edges)[, c("from", "to")]
#> from to
#> 1 subject:P1 subject:P2
#> 2 subject:P2 subject:P3
#> 3 subject:P5 subject:P6
PLINK’s --make-rel square instead writes a square GRM. Pass the matrix directly, with the subject ids as its dimnames; it is expanded to its upper-triangle pairs before thresholding, so the diagonal never becomes a self-edge:
grm <- matrix(
c(0.50, 0.25, 0.02, 0.00,
0.25, 0.50, 0.20, 0.00,
0.02, 0.20, 0.50, 0.00,
0.00, 0.00, 0.00, 0.50),
nrow = 4, byrow = TRUE,
dimnames = list(paste0("P", 1:4), paste0("P", 1:4))
)
as.data.frame(relatedness_edges_from_kinship(grm, threshold = 0.1))[
, c("from", "to")
]
#> from to
#> 1 subject:P1 subject:P2
#> 2 subject:P2 subject:P3
Whichever shape you start from, the value that passed the threshold is kept on the edge, so you can always see why a pair was linked rather than trusting the grouping blind:
kept <- as.data.frame(query_edge_type(g_rel, "subject_related_to"))
data.frame(
from = kept$from,
to = kept$to,
kinship = vapply(kept$attrs, function(a) a$kinship, numeric(1))
)
#> from to kinship
#> 1 subject:P1 subject:P2 0.25
#> 2 subject:P2 subject:P3 0.20
#> 3 subject:P5 subject:P6 0.30
The grouping is a transitive closure over the subject_related_to edges. The network below draws those edges between subjects, coloured by the relatedness group each subject (and therefore its samples) lands in: the P1–P2–P3 chain becomes one group even though P1 and P3 were never linked directly, P5–P6 form a second, and the unrelated P4 stands alone.
subject_group <- setNames(rel_groups[meta_r$sample_id], meta_r$subject_id)
kept_pairs <- kin[kin$kinship >= 0.1, c("id1", "id2")]
rel_net <- igraph::graph_from_data_frame(
kept_pairs, directed = FALSE,
vertices = data.frame(name = meta_r$subject_id)
)
palette_rel <- c("#4C78A8", "#F58518", "#54A24B", "#B279A2")
set.seed(1)
plot(rel_net,
vertex.color = palette_rel[as.integer(factor(subject_group[igraph::V(rel_net)$name]))],
vertex.size = 34,
vertex.label.color = "white",
vertex.label.font = 2,
edge.color = "grey60",
edge.width = 2,
main = "Relatedness clusters (kinship >= 0.1)")
The threshold is the key knob, and it belongs to the edge-building step, not the grouping. Raising it drops weaker links: at 0.22 the P2–P3 pair (kinship 0.20) no longer qualifies, so that chain breaks and P3 splits into its own group, while the stronger P5–P6 pair is untouched:
rel_strict <- relatedness_edges_from_kinship(kin, threshold = 0.22)
g_rel_strict <- build_dependency_graph(list(samples, subjects), list(belongs, rel_strict))
grouping_vector(derive_split_constraints(g_rel_strict, mode = "relatedness"))
#> S1 S2 S3
#> "relatedness:component_1" "relatedness:component_1" "relatedness:component_2"
#> S4 S5 S6
#> "relatedness:component_3" "relatedness:component_4" "relatedness:component_4"
Push the threshold further and the relation stops doing much work: samples end up alone in their own groups, which protects nothing. That case is detected — once more than half the samples sit in a group of one, the pairwise deriver records it in metadata$warnings, so a cut that quietly switched the relation off does not pass unnoticed:
sparse <- derive_split_constraints(
build_dependency_graph(
list(samples, subjects),
list(belongs, relatedness_edges_from_kinship(kin, threshold = 0.28))
),
mode = "relatedness"
)
sparse$metadata$n_groups
#> [1] 5
sparse$metadata$warnings
#> [1] "Most relatedness groups are singletons; pairwise coverage may be sparse or the threshold may be too strict."
That is the sparse end of the same trade-off the collapsing composite showed at the dense end: too permissive a cut merges the cohort into one unusable group, too strict a cut dissolves it into singletons that protect nothing. The threshold is where you choose between them.
Pairwise relation: spatial proximity
Spatial proximity works the same way over sample coordinates — for example spot locations from spatial transcriptomics, positions on a tissue slide, or geographic site coordinates. spatial_edges_from_coords() connects samples within a radius (Euclidean distance over the coordinate columns), and mode = "spatial" groups the resulting connected components. The distance is computed over as many coordinate columns as you give it, so a z column for a tissue volume works exactly like a plain x/y slide.
# Two spatial clusters. Cluster 1 (S1-S3) is a chain: neighbouring pairs are
# within the radius, but the endpoints are not.
coords <- data.frame(
sample_id = paste0("S", 1:6),
x = c(0, 1, 2, 6.0, 6.9, 6.2),
y = c(0, 1, 0, 6.0, 6.6, 5.3),
stringsAsFactors = FALSE
)
adj_edges <- spatial_edges_from_coords(coords, radius = 1.5)
meta_s <- data.frame(
sample_id = paste0("S", 1:6),
subject_id = paste0("P", 1:6),
stringsAsFactors = FALSE
)
samples_s <- create_nodes(meta_s, "Sample", "sample_id")
subjects_s <- create_nodes(meta_s, "Subject", "subject_id")
belongs_s <- create_edges(meta_s, "sample_id", "subject_id",
"Sample", "Subject", "sample_belongs_to_subject")
g_sp <- build_dependency_graph(list(samples_s, subjects_s), list(belongs_s, adj_edges))
sp_groups <- grouping_vector(derive_split_constraints(g_sp, mode = "spatial"))
sp_groups
#> S1 S2 S3
#> "spatial:component_1" "spatial:component_1" "spatial:component_1"
#> S4 S5 S6
#> "spatial:component_2" "spatial:component_2" "spatial:component_2"
Plotting the coordinates, drawing the within-radius adjacency edges in grey, and colouring points by the derived group makes the transitive closure concrete: S1–S2 and S2–S3 are each within the 1.5 radius, so all three share a group even though S1 and S3 are 2 units apart and were never linked directly. Every sample in the second cluster is likewise reachable from the others, while the two clusters are far enough apart to stay separate:
sp_grp <- factor(sp_groups[coords$sample_id])
row_of <- setNames(seq_len(nrow(coords)), coords$sample_id)
from_i <- row_of[sub("^sample:", "", adj_edges$data$from)]
to_i <- row_of[sub("^sample:", "", adj_edges$data$to)]
palette_sp <- c("#4C78A8", "#F58518")
plot(coords$x, coords$y, type = "n", asp = 1, xlab = "x", ylab = "y",
main = "Spatial groups (radius = 1.5)")
segments(coords$x[from_i], coords$y[from_i],
coords$x[to_i], coords$y[to_i], col = "grey60", lwd = 2)
points(coords$x, coords$y, pch = 19, cex = 3.5, col = palette_sp[as.integer(sp_grp)])
text(coords$x, coords$y, labels = coords$sample_id, col = "white", cex = 0.8, font = 2)
legend("topleft", legend = levels(sp_grp), pch = 19,
col = palette_sp[seq_along(levels(sp_grp))], title = "Spatial group", bty = "n")
One thing to watch in a real coordinate table. When coord_cols is not given, every numeric column except the id is treated as a coordinate — so a QC metric or a slide number sitting in the same frame silently joins the distance calculation and can switch the relation off entirely:
coords_qc <- cbind(coords, reads_millions = c(31, 44, 12, 38, 27, 51))
# `reads_millions` is picked up as a third dimension: nothing is within radius.
nrow(as.data.frame(spatial_edges_from_coords(coords_qc, radius = 1.5)))
#> [1] 0
# Naming the coordinate columns restores the five adjacency edges.
nrow(as.data.frame(
spatial_edges_from_coords(coords_qc, radius = 1.5, coord_cols = c("x", "y"))
))
#> [1] 5
Pass coord_cols whenever the frame carries anything numeric that is not a coordinate.
Deriving on a subset is leakage-safe
Real splits are derived on a subset of samples — the training rows, say. For pairwise (and composite) modes this raises a subtle question: if a sample that bridges two others is left out of the subset, could those two still inherit a shared group from the full graph? They do not. When you pass samples =, grouping is recomputed within that subset, so structure that exists only through an excluded sample never leaks across the split.
The spatial chain makes this visible. S1 and S3 shared a group only because S2 bridged them; ask for S1 and S3 alone, and they correctly fall into separate groups:
grouping_vector(
derive_split_constraints(g_sp, mode = "spatial", samples = c("S1", "S3"))
)
#> S1 S3
#> "spatial:component_1" "spatial:component_2"
A subset is sometimes better expressed as a graph in its own right, for example when you want to validate it, plot it, or hand it around. subset_graph() does that, and it uses the same rule, so the grouping agrees with samples =:
g_two <- subset_graph(g_sp, samples = c("S1", "S3"), graph_name = "S1 + S3")
grouping_vector(derive_split_constraints(g_two, mode = "spatial"))
#> S1 S3
#> "spatial:component_1" "spatial:component_2"
Combining a pairwise relation with a direct one
Pairwise relations are not a separate world: either of them can be listed in a composite via alongside the direct relations, and all of them feed the same connected-component search. That matters when two different mechanisms each link a different pair of samples, and severing only one would leave the other open.
Here relatedness links two subjects that share no batch, while batch links two samples from unrelated subjects. Neither relation alone separates the cohort correctly; the composite does.
The graphs earlier on this page were assembled node set by node set to keep every piece visible, but a pairwise relation does not require that. Build the direct structure from the metadata frame as usual, then attach the thresholded edges to it with add_edges():
meta_c <- data.frame(
sample_id = paste0("S", 1:4),
subject_id = paste0("P", 1:4),
batch_id = c("B1", "B1", "B2", "B3"),
stringsAsFactors = FALSE
)
kin_c <- data.frame(id1 = "P2", id2 = "P3", kinship = 0.25, stringsAsFactors = FALSE)
g_mixed <- add_edges(
graph_from_metadata(meta_c, graph_name = "mixed-structure"),
relatedness_edges_from_kinship(kin_c, threshold = 0.1)
)
# batch alone: {S1,S2} {S3} {S4} relatedness alone: {S2,S3}
grouping_vector(derive_split_constraints(g_mixed, "batch"))
#> S1 S2 S3 S4
#> "batch:B1" "batch:B1" "batch:B2" "batch:B3"
grouping_vector(derive_split_constraints(g_mixed, "relatedness"))
#> S1 S2 S3
#> "relatedness:component_1" "relatedness:component_2" "relatedness:component_2"
#> S4
#> "relatedness:component_3"
# combined, the chain S1-S2 (batch) - S3 (kinship) becomes one group:
mixed <- derive_split_constraints(g_mixed, "composite", via = c("batch", "relatedness"))
grouping_vector(mixed)
#> S1 S2 S3 S4
#> "component_1" "component_1" "component_1" "component_2"
mixed$metadata$via
#> [1] "Batch" "relatedness"
Note what via accepts and where. derive_split_constraints() takes either a lower-case mode or a node type, and pairwise modes are valid there — that is what makes the call above possible. The graph-level views take node types only: detect_dependency_components() and plot(focus = "sample_projection") project samples through Subject, Batch and friends, and reject "relatedness", because a pairwise relation is an edge set rather than a node type to route through. Use the derivation to combine them and grouping_vector() to read the result.
Under strategy = "rule_based" a pairwise source behaves as a fallback instead: a sample takes the first relation in priority that gives it a non-singleton group, so relations never chain.
Whichever route you take, the threshold that produced the edges is recorded on the graph, per relation, and is serialised with it — so a reader of the file can see which cut created the grouping:
g_mixed$metadata$edge_sources$subject_related_to[c("threshold", "metric")]
#> $threshold
#> [1] 0.1
#>
#> $metric
#> [1] "kinship"
A spec derived from a single pairwise mode also carries that threshold as a scalar in its own metadata. A composite spec does not, because several relations (each with its own cut) may have contributed; the graph’s edge_sources is the complete record in that case:
as_split_spec(derive_split_constraints(g_mixed, "relatedness"),
graph = g_mixed)$metadata[c("threshold", "threshold_metric")]
#> $threshold
#> [1] 0.1
#>
#> $threshold_metric
#> [1] "kinship"