diff --git a/DESCRIPTION b/DESCRIPTION index 88cab66..3535f31 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,7 +1,7 @@ Package: visPedigree Type: Package Title: Tidying, Analysis, and Fast Visualization of Animal and Plant Pedigrees -Version: 1.10.0 +Version: 1.10.1 Authors@R: person("Sheng","Luan", email="luansheng@gmail.com", role = c("aut","cre")) Description: Provides tools for the analysis and visualization of animal and plant pedigrees. Analytical methods include equivalent complete generations, diff --git a/NEWS.md b/NEWS.md index d943bfb..641a430 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,11 @@ +# Changes in version 1.10.1 + +## Bug fixes +1. **`pedexport(software = "wombat")` now uses the integer layout**: The 1.10.0 character-ID WOMBAT format was rejected by the real WOMBAT program — the WOMBAT manual (section 6.3) requires three integer variables with codes in 0..2147483647, offspring codes numerically larger than either parent, and unknown parents coded `0`. The `wombat` format now shares the integer layout with `blupf90` (including the `xref` attribute and `.xref` mapping file), satisfying these requirements by construction. Non-zero custom `missing` symbols are rejected for `wombat`. Verified end-to-end against the WOMBAT 26-05-2025 binary. + +## New features +1. **`pedexport(software = "hiblup")`**: New format for HIBLUP's `--pedigree` file: character columns `animal`/`sire`/`dam`, no header by default, missing parents coded `"0"`. Verified end-to-end against the HIBLUP binary (estimated additive variance and EBVs identical to ASReml/WOMBAT on the same pedigree). + # Changes in version 1.10.0 released on 8 Aug 2026 ## New features 1. **`pedexport()` for breeding software formats**: New function `pedexport()` converts a `tidyped` pedigree into the file format of common animal and plant breeding programs - BLUPF90, ASReml, Echidna, WOMBAT, MTDFREML, DMU - plus a generic `numeric` layout and an in-memory `sommer` format. Character formats (`asreml`, `echidna`, `wombat`, `sommer`) keep the original character IDs; numeric formats (`blupf90`, `mtdfreml`, `dmu`, `numeric`) renumber the pedigree and carry an `xref` mapping back to the original IDs (attribute on the returned table, and a `.xref` file when `file` is given), mirroring RENUMF90's `_XrefID` file. Defaults follow each program's conventions for header, separator, and missing-parent symbol, and every format sorts rows so parents precede offspring. diff --git a/R/pedexport.R b/R/pedexport.R index 08c2e5a..ec4b8b1 100644 --- a/R/pedexport.R +++ b/R/pedexport.R @@ -21,15 +21,23 @@ #' parents encoded as \code{"0"}. Rows are sorted by generation #' (\code{Gen}) so founders appear first. Compatible with ASReml-R #' and ASReml-SA (character IDs require the \code{!ALPHA} qualifier; -#' the header line requires the \code{!SKIP 1} qualifier).} +#' the header line requires the \code{!SKIP 1} qualifier, placed +#' immediately after the pedigree file name).} #' \item{\code{"echidna"}}{Identical character layout and defaults to #' \code{"asreml"}: columns \code{animal}, \code{sire}, and \code{dam}, #' a header by default, missing parents encoded as \code{"0"}, and rows #' sorted by generation.} -#' \item{\code{"wombat"}}{Three character columns (\code{animal}, -#' \code{sire}, \code{dam}), no header, missing parents encoded as -#' \code{"0"}. Wombat accepts alphanumeric IDs and recodes them -#' internally. Rows are sorted so parents precede offspring.} +#' \item{\code{"hiblup"}}{Three character columns (\code{animal}, +#' \code{sire}, \code{dam}), no header by default, missing parents +#' encoded as \code{"0"}, and rows sorted by generation. Compatible +#' with HIBLUP's \code{--pedigree} file.} +#' \item{\code{"wombat"}}{Integer layout identical to \code{"blupf90"} +#' (\code{IndNum}, \code{SireNum}, \code{DamNum}), no header, missing +#' parents encoded as \code{0}. WOMBAT requires integer codes in the +#' range 0 to 2147483647, with every offspring numerically larger than +#' either parent and unknown parents coded \code{0} (WOMBAT manual +#' section 6.3); the \code{tidyped()} integer index satisfies both +#' ordering constraints by construction.} #' \item{\code{"mtdfreml"}}{Identical integer layout to \code{"blupf90"}. #' Compatible with MTDFREML and MTGSAM.} #' \item{\code{"dmu"}}{Identical integer layout to \code{"blupf90"}. @@ -51,27 +59,30 @@ #' @param sep Character scalar. Field separator used when writing to #' \code{file}. Defaults to a single space (\code{" "}), which every #' file-based format accepts. BLUPF90 requires a space; WOMBAT, MTDFREML, -#' and DMU accept spaces or TABs; ASReml and Echidna accept spaces, TABs, -#' or commas; and the generic \code{"numeric"} format accepts any +#' DMU, and HIBLUP accept spaces or TABs; ASReml and Echidna accept spaces, +#' TABs, or commas; and the generic \code{"numeric"} format accepts any #' single-byte, non-empty separator supported by #' \code{\link[data.table]{fwrite}}. ASReml/Echidna comma-delimited output #' requires a \code{.csv} file extension or the \code{!CSV} qualifier. #' @param header Logical scalar or \code{NULL}. Whether to include a column #' header line. \code{NULL} (default) uses the software-specific default: #' \code{TRUE} for \code{"asreml"}, \code{"echidna"}, and \code{"numeric"}, -#' \code{FALSE} for \code{"blupf90"}, \code{"wombat"}, \code{"mtdfreml"} -#' and \code{"dmu"}. Ignored for \code{"sommer"}. Note for ASReml and -#' Echidna: a header line is read as data unless the \code{!SKIP 1} -#' qualifier is used in the command file. +#' \code{FALSE} for \code{"blupf90"}, \code{"wombat"}, \code{"mtdfreml"}, +#' \code{"dmu"}, and \code{"hiblup"}. Ignored for \code{"sommer"}. Note +#' for ASReml and Echidna: a header line is read as data unless the +#' \code{!SKIP 1} qualifier is used in the command file (it must be placed +#' immediately after the pedigree file name — trailing \code{!SKIP} after +#' other qualifiers is silently ignored by ASReml-SA). #' @param missing Character or integer scalar. Symbol for missing parents. #' \code{NULL} (default) uses the software-specific default: \code{0L} for #' numeric formats, \code{"0"} for \code{"asreml"}, \code{"echidna"}, and -#' \code{"wombat"}. -#' Numeric formats (\code{"blupf90"}, \code{"mtdfreml"}, \code{"dmu"}, -#' \code{"numeric"}) require a single integer value; \code{"asreml"}, -#' \code{"echidna"}, and \code{"wombat"} accept a character value (numeric -#' values are converted to character). Ignored for \code{"sommer"}, which -#' always codes missing parents as \code{NA}. +#' \code{"hiblup"}. +#' Numeric formats (\code{"blupf90"}, \code{"wombat"}, \code{"mtdfreml"}, +#' \code{"dmu"}, \code{"numeric"}) require a single integer value, and for +#' \code{"wombat"} the value must be \code{0} (WOMBAT manual section 6.3); +#' \code{"asreml"}, \code{"echidna"}, and \code{"hiblup"} accept a +#' character value (numeric values are converted to character). Ignored +#' for \code{"sommer"}, which always codes missing parents as \code{NA}. #' #' @return A \code{data.table} in the target format, returned invisibly. #' Numeric formats carry an \code{xref} attribute mapping each numeric ID @@ -88,9 +99,9 @@ #' \tabular{llll}{ #' \strong{software} \tab \strong{Col 1} \tab \strong{Col 2} \tab #' \strong{Col 3} \cr -#' blupf90 / mtdfreml / dmu / numeric \tab \code{IndNum} (integer) \tab -#' \code{SireNum} (integer) \tab \code{DamNum} (integer) \cr -#' asreml / echidna / wombat \tab \code{animal} (character) \tab +#' blupf90 / wombat / mtdfreml / dmu / numeric \tab \code{IndNum} +#' (integer) \tab \code{SireNum} (integer) \tab \code{DamNum} (integer) \cr +#' asreml / echidna / hiblup \tab \code{animal} (character) \tab #' \code{sire} (character) \tab \code{dam} (character) \cr #' sommer \tab \code{ID} (character) \tab \code{Sire} (character) \tab #' \code{Dam} (character) \cr @@ -107,8 +118,11 @@ #' \strong{missing} \tab \strong{IDs} \tab \strong{notes} \cr #' blupf90 \tab no \tab spaces only \tab \code{0} \tab integer \tab #' TAB separators rejected \cr -#' wombat \tab no \tab space / TAB \tab \code{"0"} \tab character \tab -#' alphanumeric accepted, recoded internally \cr +#' wombat \tab no \tab space / TAB \tab \code{0} \tab integer \tab +#' integer codes 0..2147483647; offspring code must exceed both parent +#' codes (manual section 6.3) \cr +#' hiblup \tab no \tab space / TAB \tab \code{"0"} \tab character \tab +#' alphanumeric accepted; missing parents may also be \code{NA} \cr #' mtdfreml \tab no \tab space / TAB \tab \code{0} \tab integer \tab \cr #' dmu \tab no \tab space / TAB \tab \code{0} \tab integer \tab \cr #' asreml \tab yes \tab space / TAB / comma \tab \code{"0"} \tab @@ -127,7 +141,7 @@ #' All numeric formats sort by \code{IndNum} ascending, which guarantees that #' parent rows appear before offspring rows (parents always receive a smaller #' integer index after \code{tidyped()} topological sorting). The -#' \code{"asreml"}, \code{"echidna"}, \code{"wombat"} and \code{"sommer"} +#' \code{"asreml"}, \code{"echidna"}, \code{"hiblup"} and \code{"sommer"} #' formats sort by \code{Gen} ascending for the same reason. #' #' \strong{Optional \code{tidyped()} columns:} @@ -182,6 +196,16 @@ #' out_echidna <- pedexport(tp, software = "echidna") #' identical(out_echidna, out_asreml) #' +#' # HIBLUP: character IDs, no header by default +#' out_hiblup <- pedexport(tp, software = "hiblup") +#' head(out_hiblup) +#' +#' # WOMBAT requires integer codes (manual section 6.3): same integer +#' # layout as BLUPF90, with the xref mapping back to original IDs +#' out_wombat <- pedexport(tp, software = "wombat") +#' head(out_wombat) +#' head(attr(out_wombat, "xref")) +#' #' # Numeric formats carry the ID mapping back to character IDs #' out_dmu <- pedexport(tp, software = "dmu") #' head(attr(out_dmu, "xref")) @@ -205,8 +229,9 @@ #' @import data.table #' @export pedexport <- function(ped, - software = c("blupf90", "asreml", "echidna", "wombat", - "mtdfreml", "dmu", "numeric", "sommer"), + software = c("blupf90", "asreml", "echidna", "hiblup", + "wombat", "mtdfreml", "dmu", "numeric", + "sommer"), file = NULL, sep = " ", header = NULL, @@ -239,7 +264,7 @@ pedexport <- function(ped, # ---- 2. Software-specific defaults ---- asreml_like <- c("asreml", "echidna") - use_char <- software %in% c(asreml_like, "wombat", "sommer") + use_char <- software %in% c(asreml_like, "hiblup", "sommer") def_header <- software %in% c(asreml_like, "numeric") def_miss <- if (use_char) "0" else 0L @@ -253,22 +278,27 @@ pedexport <- function(ped, } else if (use_char) { if ((!is.character(missing) && !is.numeric(missing)) || length(missing) != 1L || is.na(missing)) { - stop("For the 'asreml', 'echidna', and 'wombat' formats, 'missing' ", + stop("For the 'asreml', 'echidna', and 'hiblup' formats, 'missing' ", "must be a single character (or numeric) value.", call. = FALSE) } missing <- as.character(missing) if (!nzchar(missing)) { - stop("For the 'asreml', 'echidna', and 'wombat' formats, 'missing' ", + stop("For the 'asreml', 'echidna', and 'hiblup' formats, 'missing' ", "must not be empty.", call. = FALSE) } } else { missing_int <- suppressWarnings(as.integer(missing)) if (!is.numeric(missing) || length(missing) != 1L || is.na(missing) || is.na(missing_int) || missing_int != missing) { - stop("For numeric formats (blupf90, mtdfreml, dmu, numeric), ", + stop("For numeric formats (blupf90, wombat, mtdfreml, dmu, numeric), ", "'missing' must be a single integer value.", call. = FALSE) } missing <- missing_int + if (software == "wombat" && missing != 0L) { + stop("WOMBAT requires unknown parents to be coded 0 (manual section ", + "6.3); 'missing' must be 0 for software = \"wombat\".", + call. = FALSE) + } } # ---- 4. Build output table ---- @@ -406,7 +436,7 @@ pedexport <- function(ped, stop("BLUPF90 pedigree files require sep = \" \".", call. = FALSE) } - whitespace_formats <- c("wombat", "mtdfreml", "dmu") + whitespace_formats <- c("wombat", "mtdfreml", "dmu", "hiblup") if (software %in% whitespace_formats && !(sep %in% c(" ", "\t"))) { stop(sprintf("%s pedigree files require a space or TAB separator.", toupper(software)), call. = FALSE) @@ -431,7 +461,7 @@ pedexport <- function(ped, #' @noRd .validate_export_fields <- function(ped, software, missing, sep) { fields <- ped$Ind - if (software %in% c("asreml", "echidna", "wombat")) { + if (software %in% c("asreml", "echidna", "hiblup")) { fields <- c(fields, ped$Sire, ped$Dam, missing) } fields <- as.character(fields[!is.na(fields)]) diff --git a/man/pedexport.Rd b/man/pedexport.Rd index 376793d..c41ec5e 100644 --- a/man/pedexport.Rd +++ b/man/pedexport.Rd @@ -6,8 +6,8 @@ \usage{ pedexport( ped, - software = c("blupf90", "asreml", "echidna", "wombat", "mtdfreml", "dmu", "numeric", - "sommer"), + software = c("blupf90", "asreml", "echidna", "hiblup", "wombat", "mtdfreml", "dmu", + "numeric", "sommer"), file = NULL, sep = " ", header = NULL, @@ -31,15 +31,23 @@ and \code{Dam} identifier must be present in \code{Ind}.} parents encoded as \code{"0"}. Rows are sorted by generation (\code{Gen}) so founders appear first. Compatible with ASReml-R and ASReml-SA (character IDs require the \code{!ALPHA} qualifier; - the header line requires the \code{!SKIP 1} qualifier).} + the header line requires the \code{!SKIP 1} qualifier, placed + immediately after the pedigree file name).} \item{\code{"echidna"}}{Identical character layout and defaults to \code{"asreml"}: columns \code{animal}, \code{sire}, and \code{dam}, a header by default, missing parents encoded as \code{"0"}, and rows sorted by generation.} - \item{\code{"wombat"}}{Three character columns (\code{animal}, - \code{sire}, \code{dam}), no header, missing parents encoded as - \code{"0"}. Wombat accepts alphanumeric IDs and recodes them - internally. Rows are sorted so parents precede offspring.} + \item{\code{"hiblup"}}{Three character columns (\code{animal}, + \code{sire}, \code{dam}), no header by default, missing parents + encoded as \code{"0"}, and rows sorted by generation. Compatible + with HIBLUP's \code{--pedigree} file.} + \item{\code{"wombat"}}{Integer layout identical to \code{"blupf90"} + (\code{IndNum}, \code{SireNum}, \code{DamNum}), no header, missing + parents encoded as \code{0}. WOMBAT requires integer codes in the + range 0 to 2147483647, with every offspring numerically larger than + either parent and unknown parents coded \code{0} (WOMBAT manual + section 6.3); the \code{tidyped()} integer index satisfies both + ordering constraints by construction.} \item{\code{"mtdfreml"}}{Identical integer layout to \code{"blupf90"}. Compatible with MTDFREML and MTGSAM.} \item{\code{"dmu"}}{Identical integer layout to \code{"blupf90"}. @@ -63,8 +71,8 @@ which is intended for direct use as an R \code{data.table}.} \item{sep}{Character scalar. Field separator used when writing to \code{file}. Defaults to a single space (\code{" "}), which every file-based format accepts. BLUPF90 requires a space; WOMBAT, MTDFREML, -and DMU accept spaces or TABs; ASReml and Echidna accept spaces, TABs, -or commas; and the generic \code{"numeric"} format accepts any +DMU, and HIBLUP accept spaces or TABs; ASReml and Echidna accept spaces, +TABs, or commas; and the generic \code{"numeric"} format accepts any single-byte, non-empty separator supported by \code{\link[data.table]{fwrite}}. ASReml/Echidna comma-delimited output requires a \code{.csv} file extension or the \code{!CSV} qualifier.} @@ -72,20 +80,23 @@ requires a \code{.csv} file extension or the \code{!CSV} qualifier.} \item{header}{Logical scalar or \code{NULL}. Whether to include a column header line. \code{NULL} (default) uses the software-specific default: \code{TRUE} for \code{"asreml"}, \code{"echidna"}, and \code{"numeric"}, -\code{FALSE} for \code{"blupf90"}, \code{"wombat"}, \code{"mtdfreml"} -and \code{"dmu"}. Ignored for \code{"sommer"}. Note for ASReml and -Echidna: a header line is read as data unless the \code{!SKIP 1} -qualifier is used in the command file.} +\code{FALSE} for \code{"blupf90"}, \code{"wombat"}, \code{"mtdfreml"}, +\code{"dmu"}, and \code{"hiblup"}. Ignored for \code{"sommer"}. Note +for ASReml and Echidna: a header line is read as data unless the +\code{!SKIP 1} qualifier is used in the command file (it must be placed +immediately after the pedigree file name — trailing \code{!SKIP} after +other qualifiers is silently ignored by ASReml-SA).} \item{missing}{Character or integer scalar. Symbol for missing parents. \code{NULL} (default) uses the software-specific default: \code{0L} for numeric formats, \code{"0"} for \code{"asreml"}, \code{"echidna"}, and -\code{"wombat"}. -Numeric formats (\code{"blupf90"}, \code{"mtdfreml"}, \code{"dmu"}, -\code{"numeric"}) require a single integer value; \code{"asreml"}, -\code{"echidna"}, and \code{"wombat"} accept a character value (numeric -values are converted to character). Ignored for \code{"sommer"}, which -always codes missing parents as \code{NA}.} +\code{"hiblup"}. +Numeric formats (\code{"blupf90"}, \code{"wombat"}, \code{"mtdfreml"}, +\code{"dmu"}, \code{"numeric"}) require a single integer value, and for +\code{"wombat"} the value must be \code{0} (WOMBAT manual section 6.3); +\code{"asreml"}, \code{"echidna"}, and \code{"hiblup"} accept a +character value (numeric values are converted to character). Ignored +for \code{"sommer"}, which always codes missing parents as \code{NA}.} } \value{ A \code{data.table} in the target format, returned invisibly. @@ -108,9 +119,9 @@ result as a \code{data.table} or writes it to a plain-text file. \tabular{llll}{ \strong{software} \tab \strong{Col 1} \tab \strong{Col 2} \tab \strong{Col 3} \cr - blupf90 / mtdfreml / dmu / numeric \tab \code{IndNum} (integer) \tab - \code{SireNum} (integer) \tab \code{DamNum} (integer) \cr - asreml / echidna / wombat \tab \code{animal} (character) \tab + blupf90 / wombat / mtdfreml / dmu / numeric \tab \code{IndNum} + (integer) \tab \code{SireNum} (integer) \tab \code{DamNum} (integer) \cr + asreml / echidna / hiblup \tab \code{animal} (character) \tab \code{sire} (character) \tab \code{dam} (character) \cr sommer \tab \code{ID} (character) \tab \code{Sire} (character) \tab \code{Dam} (character) \cr @@ -127,8 +138,11 @@ below summarises what each program expects, and the defaults \strong{missing} \tab \strong{IDs} \tab \strong{notes} \cr blupf90 \tab no \tab spaces only \tab \code{0} \tab integer \tab TAB separators rejected \cr - wombat \tab no \tab space / TAB \tab \code{"0"} \tab character \tab - alphanumeric accepted, recoded internally \cr + wombat \tab no \tab space / TAB \tab \code{0} \tab integer \tab + integer codes 0..2147483647; offspring code must exceed both parent + codes (manual section 6.3) \cr + hiblup \tab no \tab space / TAB \tab \code{"0"} \tab character \tab + alphanumeric accepted; missing parents may also be \code{NA} \cr mtdfreml \tab no \tab space / TAB \tab \code{0} \tab integer \tab \cr dmu \tab no \tab space / TAB \tab \code{0} \tab integer \tab \cr asreml \tab yes \tab space / TAB / comma \tab \code{"0"} \tab @@ -147,7 +161,7 @@ below summarises what each program expects, and the defaults All numeric formats sort by \code{IndNum} ascending, which guarantees that parent rows appear before offspring rows (parents always receive a smaller integer index after \code{tidyped()} topological sorting). The -\code{"asreml"}, \code{"echidna"}, \code{"wombat"} and \code{"sommer"} +\code{"asreml"}, \code{"echidna"}, \code{"hiblup"} and \code{"sommer"} formats sort by \code{Gen} ascending for the same reason. \strong{Optional \code{tidyped()} columns:} @@ -199,6 +213,16 @@ head(out_asreml) out_echidna <- pedexport(tp, software = "echidna") identical(out_echidna, out_asreml) +# HIBLUP: character IDs, no header by default +out_hiblup <- pedexport(tp, software = "hiblup") +head(out_hiblup) + +# WOMBAT requires integer codes (manual section 6.3): same integer +# layout as BLUPF90, with the xref mapping back to original IDs +out_wombat <- pedexport(tp, software = "wombat") +head(out_wombat) +head(attr(out_wombat, "xref")) + # Numeric formats carry the ID mapping back to character IDs out_dmu <- pedexport(tp, software = "dmu") head(attr(out_dmu, "xref")) diff --git a/tests/testthat/test-pedexport-extra.R b/tests/testthat/test-pedexport-extra.R new file mode 100644 index 0000000..67772fd --- /dev/null +++ b/tests/testthat/test-pedexport-extra.R @@ -0,0 +1,212 @@ +# Supplementary systematic tests for pedexport(): format invariants across +# datasets, file round-trips, xref consistency, and custom missing symbols. +# Complements test-pedexport.R (which covers API behaviour on small examples). + +library(data.table) + +# ---- helpers ----------------------------------------------------------------- + +# Invariants that every numeric (integer-coded) export must satisfy. +expect_valid_numeric_export <- function(out, tp) { + expect_identical(names(out), c("IndNum", "SireNum", "DamNum")) + expect_true(is.integer(out$IndNum) || is.numeric(out$IndNum)) + expect_identical(out$IndNum, seq_len(nrow(tp))) + # every known parent has a strictly smaller code than the offspring + expect_true(all(out$SireNum == 0L | out$SireNum < out$IndNum)) + expect_true(all(out$DamNum == 0L | out$DamNum < out$IndNum)) + # xref maps back to the original IDs bijectively + xr <- attr(out, "xref") + expect_equal(nrow(xr), nrow(tp)) + expect_setequal(xr$Ind, tp$Ind) + expect_identical(xr$IndNum, out$IndNum) + # numeric codes identify the same parents as the character IDs + sire_char <- xr$Ind[match(out$SireNum, xr$IndNum)] + dam_char <- xr$Ind[match(out$DamNum, xr$IndNum)] + tp_ord <- tp[order(IndNum)] + expect_identical(sire_char, tp_ord$Sire) + expect_identical(dam_char, tp_ord$Dam) +} + +# Invariants for character exports: parents precede offspring, founders first, +# every parent either appears as a record or is the missing symbol. +expect_valid_char_export <- function(out, tp, miss = "0") { + expect_equal(ncol(out), 3) + ids <- out[[1]] + expect_setequal(ids, tp$Ind) + sire_pos <- match(out[[2]], ids) + dam_pos <- match(out[[3]], ids) + is_miss_s <- if (is.na(miss)) is.na(out[[2]]) else out[[2]] == miss + is_miss_d <- if (is.na(miss)) is.na(out[[3]]) else out[[3]] == miss + expect_true(all(is_miss_s | sire_pos < seq_len(nrow(out)))) + expect_true(all(is_miss_d | dam_pos < seq_len(nrow(out)))) + # structure matches the source pedigree row for row + src <- tp[match(ids, tp$Ind)] + got_s <- ifelse(is_miss_s, NA_character_, out[[2]]) + got_d <- ifelse(is_miss_d, NA_character_, out[[3]]) + expect_identical(got_s, src$Sire) + expect_identical(got_d, src$Dam) +} + +datasets <- list( + simple_ped = simple_ped, + small_ped = small_ped, + deep_ped = deep_ped, + half_founder_ped = half_founder_ped, + inbred_ped = inbred_ped +) + +# ---- 1. invariants across datasets and formats ------------------------------- + +test_that("numeric formats satisfy integer-coding invariants on all datasets", { + for (nm in names(datasets)) { + tp <- tidyped(datasets[[nm]]) + for (sw in c("blupf90", "wombat", "mtdfreml", "dmu", "numeric")) { + expect_valid_numeric_export(pedexport(tp, software = sw), tp) + } + } +}) + +test_that("character formats satisfy ordering invariants on all datasets", { + for (nm in names(datasets)) { + tp <- tidyped(datasets[[nm]]) + for (sw in c("asreml", "echidna", "hiblup")) { + expect_valid_char_export(pedexport(tp, software = sw), tp, miss = "0") + } + expect_valid_char_export(pedexport(tp, software = "sommer"), tp, miss = NA) + } +}) + +test_that("invariants hold for traced, bottom-aligned and inbreed pedigrees", { + tp1 <- tidyped(simple_ped, cand = c("J5X804", "J3Y620"), trace = "all") + tp2 <- tidyped(simple_ped, genmethod = "bottom") + tp3 <- tidyped(inbred_ped, inbreed = TRUE) + for (tp in list(tp1, tp2, tp3)) { + expect_valid_numeric_export(pedexport(tp, "blupf90"), tp) + expect_valid_char_export(pedexport(tp, "asreml"), tp, miss = "0") + } +}) + +# ---- 2. file round-trips ------------------------------------------------------ + +test_that("numeric file round-trips through fread and tidyped identically", { + tp <- tidyped(simple_ped) + f <- tempfile(fileext = ".txt") + on.exit(unlink(c(f, paste0(f, ".xref")))) + pedexport(tp, software = "blupf90", file = f) + + re <- fread(f, header = FALSE) + expect_equal(ncol(re), 3L) + expect_true(all(re$V1 == seq_len(nrow(tp)))) + tp2 <- tidyped(re) + orig <- tp[order(IndNum)] + expect_equal(nrow(tp2), nrow(tp)) + expect_identical(ifelse(is.na(tp2$Sire), 0L, as.integer(tp2$Sire)), + orig$SireNum) + expect_identical(ifelse(is.na(tp2$Dam), 0L, as.integer(tp2$Dam)), + orig$DamNum) + + # .xref file maps numeric codes back to original IDs + xr <- fread(paste0(f, ".xref"), header = FALSE) + expect_identical(xr$V1, seq_len(nrow(tp))) + expect_setequal(xr$V2, tp$Ind) +}) + +test_that("character file round-trips and preserves the pedigree structure", { + tp <- tidyped(simple_ped) + f <- tempfile(fileext = ".txt") + on.exit(unlink(f)) + pedexport(tp, software = "hiblup", file = f) # no header, "0" = missing + + re <- fread(f, header = FALSE) + expect_equal(ncol(re), 3L) + tp2 <- tidyped(re) + key <- c("Ind", "Sire", "Dam", "Gen") + expect_identical(tp2[order(Ind), ..key], tp[order(Ind), ..key]) +}) + +test_that("asreml file with header needs (and works with) an explicit skip", { + tp <- tidyped(simple_ped) + f <- tempfile(fileext = ".txt") + on.exit(unlink(f)) + pedexport(tp, software = "asreml", file = f) # header written by default + + lines <- readLines(f) + expect_match(lines[1], "animal\\s+sire\\s+dam") + re <- fread(f, skip = 1, header = FALSE) + expect_equal(nrow(re), nrow(tp)) + tp2 <- tidyped(re) + expect_setequal(tp2$Ind, tp$Ind) +}) + +# ---- 3. custom missing symbols ------------------------------------------------ + +test_that("numeric formats honour a custom integer missing symbol", { + tp <- tidyped(simple_ped) + out <- pedexport(tp, software = "dmu", missing = 9L) + founders <- tp[is.na(Sire) & is.na(Dam), IndNum] + expect_true(all(out$SireNum[founders] == 9L)) + expect_true(all(out$DamNum[founders] == 9L)) + expect_true(all(out$SireNum[-founders] != 9L | out$SireNum[-founders] == 9L & + out$SireNum[-founders] < out$IndNum[-founders])) +}) + +test_that("character formats honour a custom character missing symbol", { + tp <- tidyped(simple_ped) + out <- pedexport(tp, software = "asreml", missing = "NA") + founders <- tp[is.na(Sire) & is.na(Dam), Ind] + got <- out[match(founders, out[[1]])] + expect_true(all(got[[2]] == "NA" & got[[3]] == "NA")) + expect_valid_char_export(out, tp, miss = "NA") +}) + +# ---- 4. reconstruction paths --------------------------------------------------- + +test_that("addnum/addgen FALSE reconstruction stays valid (order may differ)", { + full <- tidyped(simple_ped) + bare <- tidyped(simple_ped, addnum = FALSE, addgen = FALSE) + + # reconstruction produces valid exports, but the row order (and hence the + # numeric coding) follows the topological row order of the input, which is + # NOT the (Gen, Ind) order used when the columns are present + num_bare <- pedexport(bare, "blupf90") + expect_identical(num_bare$IndNum, seq_len(nrow(bare))) + expect_true(all(num_bare$SireNum == 0L | num_bare$SireNum < num_bare$IndNum)) + expect_true(all(num_bare$DamNum == 0L | num_bare$DamNum < num_bare$IndNum)) + xr <- attr(num_bare, "xref") + src <- bare[match(xr$Ind, bare$Ind)] + expect_identical(xr$Ind[match(num_bare$SireNum, xr$IndNum)] , src$Sire) + expect_identical(xr$Ind[match(num_bare$DamNum, xr$IndNum)], src$Dam) + + chr_bare <- pedexport(bare, "asreml") + expect_valid_char_export(chr_bare, bare, miss = "0") + + # caller's object is not modified by reconstruction + expect_false(any(c("IndNum", "Gen") %in% names(bare))) +}) + +# ---- 5. scale ----------------------------------------------------------------- + +test_that("pedexport stays correct on a 50k simulated pedigree", { + set.seed(7) + n_founders <- 500 + n_gen <- 8 + per_gen <- 6000 + inds <- paste0("F", seq_len(n_founders)) + sex <- rep(c("male", "female"), length.out = n_founders) + recs <- list(data.table(Ind = inds, Sire = NA_character_, Dam = NA_character_)) + males <- inds[sex == "male"]; females <- inds[sex == "female"] + for (g in seq_len(n_gen)) { + new <- paste0("G", g, "_", seq_len(per_gen)) + recs[[g + 1L]] <- data.table(Ind = new, + Sire = sample(males, per_gen, replace = TRUE), + Dam = sample(females, per_gen, replace = TRUE)) + sx <- sample(c("male", "female"), per_gen, replace = TRUE) + males <- c(males, new[sx == "male"]); females <- c(females, new[sx == "female"]) + } + tp <- tidyped(rbindlist(recs)[sample(.N)]) + + out <- pedexport(tp, "blupf90") + expect_valid_numeric_export(out, tp) + outc <- pedexport(tp, "asreml") + expect_valid_char_export(outc, tp, miss = "0") +}) diff --git a/tests/testthat/test-pedexport.R b/tests/testthat/test-pedexport.R index e656406..649ecc9 100644 --- a/tests/testthat/test-pedexport.R +++ b/tests/testthat/test-pedexport.R @@ -51,23 +51,48 @@ test_that("pedexport echidna is identical to asreml", { expect_equal(echidna, asreml) }) -test_that("pedexport wombat returns character columns with '0' missing", { +test_that("pedexport wombat uses the integer layout (WOMBAT manual 6.3)", { tp <- tidyped(make_simple_ped()) womb <- pedexport(tp, software = "wombat") - # Wombat accepts alphanumeric IDs (recoding them internally), so the - # export keeps character IDs rather than renumbering. - expect_named(womb, c("animal", "sire", "dam")) - expect_type(womb$animal, "character") - expect_type(womb$sire, "character") + # WOMBAT requires integer codes 0..2147483647, offspring code larger than + # either parent, unknown parents coded 0 (manual section 6.3). + expect_named(womb, c("IndNum", "SireNum", "DamNum")) + expect_type(womb$IndNum, "integer") + expect_identical(womb$IndNum, seq_len(nrow(tp))) + expect_true(all(womb$SireNum == 0L | womb$SireNum < womb$IndNum)) + expect_true(all(womb$DamNum == 0L | womb$DamNum < womb$IndNum)) + + founders <- tp[is.na(Sire) & is.na(Dam), IndNum] + expect_true(all(womb$SireNum[founders] == 0L)) + expect_true(all(womb$DamNum[founders] == 0L)) + + # numeric formats carry the xref mapping back to original IDs + xr <- attr(womb, "xref") + expect_equal(nrow(xr), nrow(tp)) + expect_identical(xr$Ind[xr$IndNum], tp$Ind[order(tp$IndNum)]) +}) + +test_that("pedexport hiblup returns character columns with '0' missing", { + tp <- tidyped(make_simple_ped()) + hib <- pedexport(tp, software = "hiblup") + + # HIBLUP's --pedigree file: character IDs, no header by default, "0" for + # unknown parents. + expect_named(hib, c("animal", "sire", "dam")) + expect_type(hib$animal, "character") + expect_type(hib$sire, "character") founders <- tp[is.na(Sire) & is.na(Dam), Ind] - expect_true(all(womb[animal %in% founders, sire] == "0")) - expect_true(all(womb[animal %in% founders, dam] == "0")) + expect_true(all(hib[animal %in% founders, sire] == "0")) + expect_true(all(hib[animal %in% founders, dam] == "0")) # Non-missing parents retain their character IDs - expect_equal(womb[animal == "C", sire], "A") - expect_equal(womb[animal == "C", dam], "B") + expect_equal(hib[animal == "C", sire], "A") + expect_equal(hib[animal == "C", dam], "B") + + # character format: no xref attribute + expect_null(attr(hib, "xref")) }) test_that("pedexport mtdfreml is identical to blupf90", { @@ -77,6 +102,13 @@ test_that("pedexport mtdfreml is identical to blupf90", { expect_equal(b90, mtd) }) +test_that("pedexport wombat is identical to blupf90", { + tp <- tidyped(make_simple_ped()) + b90 <- pedexport(tp, software = "blupf90") + womb <- pedexport(tp, software = "wombat") + expect_equal(b90, womb) +}) + test_that("pedexport numeric returns integer columns with header by default", { tp <- tidyped(make_simple_ped()) out <- pedexport(tp, software = "numeric") @@ -201,15 +233,34 @@ test_that("echidna file uses ASReml defaults", { expect_equal(length(lines), nrow(tp) + 1L) }) -test_that("wombat file has no header by default", { +test_that("wombat file has no header by default and writes an .xref file", { tp <- tidyped(make_simple_ped()) tmp <- tempfile(fileext = ".txt") - on.exit(unlink(tmp), add = TRUE) + on.exit(unlink(c(tmp, paste0(tmp, ".xref"))), add = TRUE) pedexport(tp, software = "wombat", file = tmp) lines <- readLines(tmp) - # Wombat reads the pedigree in free format; a header line would be + # WOMBAT reads the pedigree as three integer columns; a header line would + # be parsed as data, so none is written by default. + expect_equal(length(lines), nrow(tp)) + expect_true(all(grepl("^[0-9]+ [0-9]+ [0-9]+$", lines))) + + xr_file <- paste0(tmp, ".xref") + expect_true(file.exists(xr_file)) + xr <- read.table(xr_file, header = FALSE, stringsAsFactors = FALSE) + expect_equal(nrow(xr), nrow(tp)) +}) + +test_that("hiblup file has no header by default", { + tp <- tidyped(make_simple_ped()) + tmp <- tempfile(fileext = ".txt") + on.exit(unlink(tmp), add = TRUE) + + pedexport(tp, software = "hiblup", file = tmp) + + lines <- readLines(tmp) + # HIBLUP reads the pedigree in free format; a header line would be # parsed as data, so none is written by default. expect_equal(length(lines), nrow(tp)) expect_false(lines[[1]] == "animal sire dam") @@ -257,6 +308,8 @@ test_that("pedexport validates separators by software format", { regexp = "require") expect_error(pedexport(tp, software = "wombat", sep = ","), regexp = "space or TAB") + expect_error(pedexport(tp, software = "hiblup", sep = ","), + regexp = "space or TAB") expect_error(pedexport(tp, software = "asreml", sep = "|"), regexp = "space, TAB, or comma") expect_error(pedexport(tp, software = "echidna", sep = "|"), @@ -374,6 +427,17 @@ test_that("numeric formats reject non-integer missing symbols", { expect_true(all(out[IndNum %in% founders, DamNum] == -99L)) }) +test_that("wombat requires missing parents coded 0 (manual 6.3)", { + tp <- tidyped(make_simple_ped()) + expect_error(pedexport(tp, software = "wombat", missing = -99L), + regexp = "coded 0") + expect_error(pedexport(tp, software = "wombat", missing = 9L), + regexp = "coded 0") + # the default and an explicit 0 are accepted + expect_no_error(pedexport(tp, software = "wombat")) + expect_no_error(pedexport(tp, software = "wombat", missing = 0L)) +}) + # --------------------------------------------------------------------------- # 9. dmu and sommer formats; numeric ID mapping (xref) # --------------------------------------------------------------------------- @@ -506,7 +570,7 @@ test_that("character missing symbols must be non-empty and file-safe", { test_that("numeric exports carry an xref mapping attribute", { tp <- tidyped(make_simple_ped()) - for (sw in c("blupf90", "mtdfreml", "dmu", "numeric")) { + for (sw in c("blupf90", "wombat", "mtdfreml", "dmu", "numeric")) { out <- pedexport(tp, software = sw) xr <- attr(out, "xref") expect_s3_class(xr, "data.table") @@ -520,7 +584,7 @@ test_that("character formats carry no xref attribute", { tp <- tidyped(make_simple_ped()) expect_null(attr(pedexport(tp, software = "asreml"), "xref")) expect_null(attr(pedexport(tp, software = "echidna"), "xref")) - expect_null(attr(pedexport(tp, software = "wombat"), "xref")) + expect_null(attr(pedexport(tp, software = "hiblup"), "xref")) expect_null(attr(pedexport(tp, software = "sommer"), "xref")) })