--- title: "5. Calculation and visualization of relationship matrix" author: "Sheng Luan" date: "`r Sys.Date()`" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{5. Calculation and visualization of relationship matrix} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 6.5, fig.height = 6.5, dpi = 300, out.width = "100%" ) library(visPedigree) library(Matrix) ``` 1. [Calculating Relationship Matrices with pedmat()](#1) 1.1 [Supported Methods](#1-1) 1.2 [Basic Usage](#1-2) 1.3 [Sparse Matrix Representation](#1-3) 2. [Matrix-Free Products with pedprod()](#2) 2.1 [When to Use pedprod vs pedmat](#2-1) 2.2 [Basic Usage — $Ax$ (vector)](#2-2) 2.3 [Products with $A^{-1}$](#2-3) 2.4 [Multiple Contribution Schemes — $AX$](#2-4) 2.5 [Practical Use Cases](#2-5) 2.6 [Performance and Scalability](#2-6) 3. [Inspecting the Matrix](#3) 3.1 [Summary Statistics](#3-1) 3.2 [Querying Specific Relationships](#3-2) 4. [Compact Mode for Large Pedigrees](#4) 4.1 [Using compact = TRUE](#4-1) 4.2 [Expanding and Querying Compacted Matrices](#4-2) 4.3 [When to Use Compact Mode](#4-3) 5. [Visualizing Relationship Matrices with vismat()](#5) 5.1 [Relationship Heatmaps](#5-1) 5.2 [Inbreeding and Kinship Histograms](#5-2) 6. [Performance Considerations](#6) Relationship matrices are fundamental tools in quantitative genetics and animal breeding. They quantify the genetic similarity between individuals due to shared ancestry, which is essential for estimating breeding values (BLUP) and managing genetic diversity. The `visPedigree` package provides efficient tools for calculating various relationship matrices and visualizing them through heatmaps and histograms. ## 1. Calculating Relationship Matrices with `pedmat()` {#1} The `pedmat()` function is the primary tool for calculating relationship matrices. It supports both additive and dominance relationship matrices, as well as their inverses. The inverse additive relationship matrix (`Ainv`) follows Henderson's rules (Henderson, 1976), while inbreeding coefficients use the same optimized recursive engine as `inbreed()`, based on algorithms for large populations described by Meuwissen & Luo (1992) and Sargolzaei et al. (2005). ### 1.1 Supported Methods {#1-1} The `method` parameter in `pedmat()` determines the type of matrix to calculate: - **"A"**: Additive relationship matrix (Numerator Relationship Matrix). - **"Ainv"**: Inverse of the additive relationship matrix. - **"D"**: Dominance relationship matrix. - **"Dinv"**: Inverse of the dominance relationship matrix. - **"AA"**: Additive-by-additive (epistatic) relationship matrix. - **"AAinv"**: Inverse of the epistatic relationship matrix. - **"f"**: Inbreeding coefficients vector (uses the same optimized engine as `tidyped(..., inbreed = TRUE)`). ### 1.2 Basic Usage {#1-2} Most calculations require a pedigree tidied by `tidyped()`. ```{r basic_calc} # Load example pedigree and tidy it data(small_ped) tped <- tidyped(small_ped) # Calculate Additive Relationship Matrix (A) mat_A <- pedmat(tped, method = "A") # Calculate Dominance Relationship Matrix (D) mat_D <- pedmat(tped, method = "D") # Calculate inbreeding coefficients (f) vec_f <- pedmat(tped, method = "f") ``` ### 1.3 Sparse Matrix Representation {#1-3} By default, `pedmat()` returns a sparse matrix (class `dsCMatrix` from the `Matrix` package) for relationship matrices. This is highly memory-efficient for large pedigrees where many individuals are unrelated. ```{r sparse_check} class(mat_A) ``` ## 2. Matrix-Free Products with `pedprod()` {#2} For small and moderate pedigrees, matrices returned by `pedmat()` are ordinary matrix-like R objects and can be used directly with R's `%*%`. For large pedigrees, however, forming the dense additive relationship matrix can be the dominant memory cost. `pedprod()` instead computes $Ax$, $AX$, $A^{-1}x$, or $A^{-1}X$ directly from the pedigree — without ever constructing the $n \times n$ relationship matrix. The implementation uses the pedigree factorization $A = TDT'$ (Colleau, 2002), where $T$ is the transmission matrix and $D$ is the diagonal matrix of Mendelian sampling variances. The product is computed through backward and forward pedigree traversals in $O(np)$ time for an $n \times p$ right-hand side, with $O(np)$ working memory — a fraction of the $O(n^2)$ storage that a dense $A$ would require. ### 2.1 When to Use `pedprod` vs `pedmat` {#2-1} | Use `pedprod()` when ... | Use `pedmat()` when ... | |:---|:---| | You only need matrix-vector products ($Ax$, $A^{-1}x$) | You need individual relationship coefficients | | The pedigree has >10 000 individuals | The pedigree is small enough that $A$ fits in memory | | You are iterating (e.g., BLUP, OCS, MCMC) | You need the matrix for visualization or inspection | | Memory is the bottleneck | You need $D$, $AA$, or dominance matrices | | You are working with $A^{-1}$ products for mixed models | You need to query specific pairwise relationships | The two functions are complementary: `pedmat()` gives you the full matrix when you need to examine it; `pedprod()` applies the matrix implicitly when you only need its action on vectors. ### 2.2 Basic Usage — $A x$ (vector) {#2-2} The most common use is computing the additive relationship of every individual to a weighted group — for example, a set of selection candidates. ```{r pedprod_Ax} # Load example pedigree and tidy it data(small_ped) tped <- tidyped(small_ped) # Equal contributions from two candidates; other individuals have zero weight weights <- c(Z1 = 0.5, Z2 = 0.5) # Relationship of every pedigree individual to the weighted candidate group Ax <- pedprod(tped, weights) head(Ax) # Average additive relationship c' A c and average coancestry mean_relationship <- sum(weights * Ax[names(weights)]) mean_coancestry <- mean_relationship / 2 c(mean_relationship = mean_relationship, mean_coancestry = mean_coancestry) ``` **Named vs unnamed vectors.** When `x` is named, individuals not listed automatically receive value zero — you do not need to pad the vector yourself. An unnamed vector must have exactly one entry per pedigree individual in `IndNum` order: ```{r pedprod_named_unnamed} # Named: only listed individuals get non-zero values named_x <- c(Z1 = 0.3, Z2 = 0.4, A = 0.3) pedprod(tped, named_x)[1:6] # B, C, etc. are automatically zero # Unnamed: must match the pedigree size unnamed_x <- rep(1, nrow(tped)) length(pedprod(tped, unnamed_x)) ``` Input validation catches common mistakes: duplicate or unknown IDs, missing values, and non-finite entries are all rejected with informative messages. ### 2.3 Products with $A^{-1}$ {#2-3} The inverse additive relationship matrix $A^{-1}$ is fundamental in mixed model equations (Henderson's BLUP). `pedprod()` computes $A^{-1}x$ and $A^{-1}X$ efficiently using Henderson's rules without forming the $n \times n$ matrix: ```{r pedprod_Ainv} # Ainv * vector x <- rnorm(nrow(tped)) Ainv_x <- pedprod(tped, x, method = "Ainv") head(Ainv_x) # Verify against explicit computation (small pedigree only) A <- pedmat(tped, method = "A", sparse = FALSE) Ainv <- pedmat(tped, method = "Ainv", sparse = FALSE) all.equal( unname(Ainv_x), unname(drop(Ainv %*% x)), tolerance = 1e-12 ) ``` For large pedigrees where explicit inversion is impossible, `pedprod(tped, x, method = "Ainv")` remains computationally feasible. The numerical equivalence with the full matrix can be verified on a small pedigree: ```{r pedprod_Ainv_identity} # Ainv * (A * x) should recover x (to machine precision) A_x <- pedprod(tped, x, method = "A") Ainv_Ax <- pedprod(tped, A_x, method = "Ainv") all.equal(unname(Ainv_Ax), unname(x), tolerance = 1e-12) ``` ### 2.4 Multiple Contribution Schemes — $A X$ {#2-4} When `x` is a matrix, `pedprod()` computes all columns in a single backward/ forward traversal of the pedigree. This is significantly more efficient than calling `pedprod()` repeatedly for each column. ```{r pedprod_AX} # Three ways to weight the SAME two candidates (Z1 and Z2 are full sibs) schemes <- cbind( Equal = c(Z1 = 0.5, Z2 = 0.5), Z1_only = c(Z1 = 1.0, Z2 = 0.0), Weighted = c(Z1 = 0.7, Z2 = 0.3) ) AX <- pedprod(tped, schemes) # The schemes diverge at the candidates themselves (and their descendants); # a shared ancestor is equally related to any weighting of two full sibs. AX[c("Z1", "Z2"), ] # Expected average relationship within each contributing group, 0.5 * c'A c. # Under random mating this equals the mean inbreeding of the resulting progeny. group_coancestry <- 0.5 * colSums(schemes * AX[rownames(schemes), ]) round(group_coancestry, 5) ``` The three schemes draw on the same two full sibs, yet the genetic consequence differs sharply. Concentrating all contribution on a single animal (`Z1_only`) exposes its full, partly inbred genome and inflates the group's expected progeny inbreeding to `r round(group_coancestry["Z1_only"], 3)`, whereas spreading contribution evenly (`Equal`) drives it down to `r round(group_coancestry["Equal"], 3)`. Trading a little short-term gain for a lower average relationship is exactly the balance struck by optimal contribution selection. `pedprod()` supplies the underlying $c'A c$ evaluation for every candidate column in a single $O(np)$ traversal — compared to $O(p \cdot n^2)$ for explicit products with a precomputed $A$ — and never forms $A$ itself. ### 2.5 Practical Use Cases {#2-5} #### Optimal Contribution Selection (OCS) OCS maximizes genetic gain while constraining inbreeding. The core computation repeatedly evaluates $c'A c$ for candidate contribution vectors $c$: ```{r pedprod_ocs} # Weighted candidate contributions candidates <- setNames(c(0.25, 0.25, 0.25, 0.25), c("Z1", "Z2", "A", "B")) Ac <- pedprod(tped, candidates) # Average coancestry of the selected group: 0.5 * c' A c c_accepted <- sum(candidates * Ac[names(candidates)]) / 2 c_accepted ``` With `pedprod()`, this evaluation is fast enough to be embedded inside an iterative optimization loop over thousands of candidate configurations. #### Mixed Model Equations The mixed model equations for BLUP require products $A^{-1} u$ and $A^{-1} Z$. For an $n \times p$ design matrix $Z$, use: ```{r pedprod_blup} # Simulated breeding values as a matrix right-hand side set.seed(20260704) Z_design <- cbind( trait1 = rnorm(nrow(tped)), trait2 = rnorm(nrow(tped)) ) rownames(Z_design) <- tped$Ind # Ainv * Z in one traversal — no Ainv ever stored Ainv_Z <- pedprod(tped, Z_design, method = "Ainv") dim(Ainv_Z) Ainv_Z[1:5, ] ``` #### Founder Contribution Analysis Isolating one founder at a time turns `pedprod()` into a founder-origin decomposition: the product $A e_f$ (a unit vector on founder $f$) returns that founder's expected genetic contribution to every individual. Stacking the unit vectors into a matrix decomposes the whole population in one traversal, exposing which founders dominate the current gene pool and which have been lost: ```{r pedprod_geneflow} # One unit contribution vector per founder: column f isolates founder f founder_ids <- tped[Gen == 1, Ind] founder_design <- diag(length(founder_ids)) dimnames(founder_design) <- list(founder_ids, founder_ids) # A %*% e_f is founder f's expected genetic contribution to every individual footprint <- pedprod(tped, founder_design, method = "A") # Founder composition of one candidate: contributions sum to one round(footprint["Z1", ], 4) # Average founder contribution to the youngest generation, ranked young <- tped[Gen == max(Gen), Ind] founder_share <- sort(colMeans(footprint[young, , drop = FALSE]), decreasing = TRUE) round(founder_share, 4) ``` Reading across a row (`footprint["Z1", ]`) decomposes an individual's genome into founder origins that sum to one; reading down a column isolates a single founder's footprint across the pedigree. Here founder **N** accounts for more than half of the youngest generation, while founders **I** and **R** have already fallen to zero — their founder genomes are lost. Such imbalance lowers the effective number of founders and erodes genetic diversity faster than the census size implies. Because `pedprod()` never forms $A$, the same decomposition scales to pedigrees with millions of individuals, where it underpins the monitoring of founder representation and effective population size. ### 2.6 Performance and Scalability {#2-6} The key advantage of `pedprod()` is scalability. The table below summarizes the computational cost for a pedigree with $n$ individuals and a right-hand side with $p$ columns: | $n$ | $p$ | `pedmat("A")` memory | `pedprod()` memory | Feasible with `pedmat`? | |:---|---:|:---|:---|:---| | $10^3$ | 1 | ~8 MB | ~8 KB | Yes | | $10^4$ | 1 | ~800 MB | ~80 KB | Marginal | | $10^4$ | 10 | ~800 MB | ~800 KB | Marginal | | $10^5$ | 1 | ~80 GB | ~800 KB | **No** | | $10^5$ | 100 | ~80 GB | ~8 MB | **No** | | $10^6$ | 1 | ~8 TB | ~8 MB | **No** | For pedigrees exceeding ~25 000 individuals, `pedmat(method = "A")` will refuse to construct a dense $A$ matrix. `pedprod()` has no such limit — it continues to work as long as the pedigree and right-hand side fit in memory: ```{r pedprod_large_proof, eval=FALSE} # Pedigrees beyond the dense-A guard still work with pedprod() n <- 50000L ids <- paste0("I", seq_len(n)) raw <- data.frame( Ind = ids, Sire = c(NA_character_, ids[-n]), Dam = NA_character_, stringsAsFactors = FALSE ) tped_large <- tidyped(raw) # pedprod works; pedmat would error result <- pedprod(tped_large, setNames(1, tail(ids, 1))) length(result) # 50000 ``` **Rule of thumb:** Use `pedprod()` whenever the analysis only requires a matrix product. Construct $A$, compact matrices, or grouped summaries only when the corresponding matrix entries are themselves required. ## 3. Inspecting the Matrix {#3} ### 3.1 Summary Statistics {#3-1} Use the `summary()` method to get an overview of the calculated matrix, including size, density, and average relationship. ```{r matrix_summary} tail(summary(mat_A),10) ``` ### 3.2 Querying Specific Relationships {#3-2} Instead of manually indexing the matrix, you can use `query_relationship()` to retrieve coefficients by individual IDs. ```{r query} # Query relationship between Z1 and Z2 query_relationship(mat_A, "Z1", "Z2") # Query multiple pairs query_relationship(mat_A, c("Z1", "A"), c("Z2", "B")) ``` ## 4. Compact Mode for Large Pedigrees {#4} For large pedigrees with many full-sibling families (common in aquatic breeding populations), `pedmat()` can merge full siblings into representative nodes to save memory and time. ### 4.1 Using `compact = TRUE` {#4-1} When `compact = TRUE`, the matrix is calculated for unique representative individuals from each full-sib family. ```{r compact_calc} # Calculate compacted A matrix mat_compact <- pedmat(tped, method = "A", compact = TRUE) # The result is a 'pedmat' object containing the compacted matrix print(mat_compact[11:20,11:20]) ``` ### 4.2 Expanding and Querying Compacted Matrices {#4-2} If you need the full matrix after a compact calculation, use `expand_pedmat()`. For retrieving specific values, `query_relationship()` handles both standard and compact objects transparently. ```{r expand} # Expand to full 28x28 matrix mat_full <- expand_pedmat(mat_compact) dim(mat_full) # Query still works the same way query_relationship(mat_compact, "Z1", "Z2") ``` ### 4.3 When to Use Compact Mode {#4-3} Compact mode is highly recommended for: * **Large Pedigrees**: More than 5,000 individuals with substantial full-sibling groups. * **High-fecundity species**: Such as aquatic animals or plants, where families often have hundreds or thousands of offspring. * **Memory-limited environments**: When the full matrix exceeds available RAM. | Pedigree Size | Full-Sib Proportion | Recommended Mode | | :--- | :--- | :--- | | < 1,000 | Any | Standard | | > 5,000 | < 20% | Standard / Compact | | > 5,000 | > 20% | **Compact** | ## 5. Visualizing Relationship Matrices with `vismat()` {#5} Visualization helps in understanding population structure, detecting family clusters, and checking the distribution of genetic relationships. ### 5.1 Relationship Heatmaps {#5-1} The `"heatmap"` type (default) uses a Nature Genetics style color palette (White–Orange–Red) to display relationships. Rows and columns are reordered by hierarchical clustering (Ward.D2) by default, bringing closely related individuals into contiguous blocks — full-sibs cluster tightly because they share nearly identical relationship profiles with the rest of the population. ```{r heatmap, fig.width=6, fig.height=6} # Heatmap of the A matrix (with default clustering reorder) vismat(mat_A, labelcex = 0.5) ``` #### Compact Matrix — Direct Visualization A compact `pedmat` object can be passed directly to `vismat()`. It is automatically expanded to full dimensions before rendering. ```{r heatmap_compact, fig.width=6, fig.height=6} # Compact matrix: expanded automatically (message printed) vismat(mat_compact,labelcex=0.5) ``` #### Preserve Pedigree Order Set `reorder = FALSE` to keep the original pedigree order instead of re-sorting by clustering. ```{r heatmap_no_reorder, fig.width=6, fig.height=6} vismat(mat_A, reorder = FALSE, labelcex = 0.5) ``` #### Display a Subset of Individuals Use `ids` to focus on specific individuals. ```{r heatmap_ids, fig.width=5, fig.height=5} target_ids <- rownames(as.matrix(mat_A))[1:8] vismat(mat_A, ids = target_ids, main = "Relationship Heatmap — First 8 Individuals") ``` #### Grouping by Pedigree Column For large populations, aggregate relationships to a group-level view using the `by` parameter. The matrix is reduced to mean coefficients between groups. ```{r heatmap_group, fig.width=6, fig.height=6} # Mean relationship between generations vismat(mat_A, ped = tped, by = "Gen", main = "Mean Relationship Between Generations") ``` ```{r heatmap_family, fig.width=6, fig.height=6} # Mean relationship between full-sib families # (founders without a family assignment are excluded automatically) vismat(mat_A, ped = tped, by = "Family", main = "Mean Relationship Between Full-Sib Families") ``` ### 5.2 Inbreeding and Kinship Histograms {#5-2} The "histogram" type displays the distribution of relationship coefficients (lower triangle) or inbreeding coefficients. ```{r histogram, fig.width=6, fig.height=4} # Distribution of relationship coefficients vismat(mat_A, type = "histogram") ``` ## 6. Performance Considerations {#6} Calculation and visualization of large matrices can be resource-intensive. `vismat()` applies the following automatic optimizations: | Condition | Behavior | | :--- | :--- | | Compact + `by` | Group means are computed directly from the compact matrix (no full expansion) | | Compact, no `by`, N > 5 000 | Uses compact representative view (labels show `ID (×n)`) | | Compact, no `by`, N ≤ 5 000 | Matrix is automatically expanded via `expand_pedmat()` | | N > 2 000 | Hierarchical clustering (reorder) is automatically skipped | | N > 500 | Individual labels are automatically hidden | | N > 100 | Grid lines are automatically hidden | When a compact `pedmat` is used with `by`, `vismat()` computes the group-level mean relationship matrix algebraically from the K×K compact matrix, including a sibling off-diagonal correction. This avoids expanding to the full N×N matrix, making family-level or generation-level visualization feasible even for pedigrees with hundreds of thousands of individuals. The example below uses `big_family_size_ped` (178 431 individuals, compact to 2 626) and displays the mean additive relationship among **all** full-sib families in the latest generation — a computation that would be infeasible with full expansion. ```{r large_ped_tip, fig.width=9, fig.height=8} data(big_family_size_ped) tp_big <- tidyped(big_family_size_ped) last_gen <- max(tp_big$Gen, na.rm = TRUE) # Compute the compact A matrix for the entire pedigree mat_big_compact <- pedmat(tp_big, method = "A", compact = TRUE) # Focus on all individuals in the last generation that belong to a family ids_last_gen <- tp_big[Gen == last_gen & !is.na(Family), Ind] # vismat() aggregates directly from the compact matrix — no expansion needed vismat( mat_big_compact, ped = tp_big, ids = ids_last_gen, by = "Family", labelcex = 0.3, main = paste("Mean Relationship Between All Families in Generation", last_gen) ) ``` This family-level view reveals the genetic structure among all `r length(unique(tp_big[Gen == last_gen & !is.na(Family), Family]))` families comprising `r length(ids_last_gen)` individuals, computed in seconds from the compact matrix. --- **See Also:** - `vignette("tidy-pedigree", package = "visPedigree")` - `vignette("draw-pedigree", package = "visPedigree")` ## References - Colleau, J. J. (2002). An indirect approach to the extensive calculation of relationship coefficients. *Genetics Selection Evolution*, 34, 409-421. - Henderson, C. R. (1976). A simple method for computing the inverse of a numerator relationship matrix used in prediction of breeding values. *Biometrics*, 32(1), 69-83. - Meuwissen, T. H. E., & Luo, Z. (1992). Computing inbreeding coefficients in large populations. *Genetics Selection Evolution*, 24, 305-313. - Sargolzaei, M., Iwaisaki, H., & Colleau, J.-J. (2005). A fast algorithm for computing inbreeding coefficients in large populations. *Journal of Animal Breeding and Genetics*, 122, 325-331.