From 589f6802a6b4155e2f6db2b0962fd6ad72936aba Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Mon, 12 Oct 2026 16:46:38 -0400 Subject: [PATCH 1/5] submitted 2.2 --- CRAN-SUBMISSION | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/CRAN-SUBMISSION b/CRAN-SUBMISSION index 67aa439..6f78db6 100644 --- a/CRAN-SUBMISSION +++ b/CRAN-SUBMISSION @@ -1,3 +1,3 @@ -Version: 2.1.0 -Date: 2026-10-05 14:23:25 UTC -SHA: 4ce6c92ad91c416a6d79046d7519ce2afd5b8561 +Version: 2.2.0 +Date: 2026-10-12 20:45:10 UTC +SHA: ca36034f20aec4bc46ba0a324901f7b564588e79 From 322b42f45b85dc40f18abaa7a7ffbd7be0bb3e2d Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Mon, 26 Oct 2026 14:40:18 -0400 Subject: [PATCH 2/5] Deep check of code with AI (Claude Opus 5.5) to improve reporting of cycle and dependencies on pedigree cleaning. All output reviewed by a human. Additional improvements to tutorial for clarity and thoroughness --- DESCRIPTION | 2 +- NEWS.md | 2 + R/check_ped.R | 138 ++-- docs/BIGpopA_user_guide.Rmd | 396 +++++++--- docs/BIGpopA_user_guide.html | 1229 ++++++++++++++++++----------- man/check_ped.Rd | 5 +- tests/testthat/test-check_ped.R | 68 ++ tests/testthat/test-hard-cycles.R | 301 +++++++ 8 files changed, 1523 insertions(+), 618 deletions(-) create mode 100644 tests/testthat/test-hard-cycles.R diff --git a/DESCRIPTION b/DESCRIPTION index c55da55..a75b902 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -48,4 +48,4 @@ Suggests: withr VignetteBuilder: knitr Config/testthat/edition: 3 -Config/roxygen2/version: 8.1.0 +Config/roxygen2/version: 8.1.1 diff --git a/NEWS.md b/NEWS.md index 56ae501..8641432 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,6 @@ # BIGpopA 2.2.0 +- `check_ped()`: each cycle in `dependencies` is now reported once, as a closed path in which every id is a recorded parent of the one before it (e.g. `A -> B -> A`), instead of once per individual in the cycle. Interlocked loops are reported together, with one loop as an example. Cycle detection now runs in linear time without recursion, so large pedigrees with many shared ancestors no longer stall, and an individual recorded as its own parent is reported only by its own message. +- `check_ped()`: `dependencies` always has a character `dependency` column; with no dependencies it was a data frame without columns. - `validate_pedigree()` and `find_parentage()`: `founders_file` now also accepts a character vector of ids or a data.frame (first column), in addition to a file path. All inputs can now be supplied from memory. - `validate_pedigree()`: removed the suggested replacement parents (`best_male_candidate`, `best_male_candidate_error_pct`, `best_female_candidate`, `best_female_candidate_error_pct`) from `full_results` and the per-status tables. Use `find_parentage(method = "fill_pedigree")` on `corrected_pedigree` to fill in missing parents. - `validate_pedigree()`: new `marker_summary` element in the returned list gives, for every marker, the number of trios tested, the number with a Mendelian mismatch and the mismatch percentage, sorted worst first. It covers all trios with both parents known and genotyped, and helps spot markers (e.g. paralogous or miscalled) that cause many failures. diff --git a/R/check_ped.R b/R/check_ped.R index 5099d37..382e4d0 100644 --- a/R/check_ped.R +++ b/R/check_ped.R @@ -23,7 +23,10 @@ #' \item{conflicting_trios}{IDs with conflicting male_parent or female_parent assignments.} #' \item{inconsistent_sex_roles}{Rows where a conflicting ID appears as male_parent or female_parent.} #' \item{missing_parents}{Parent IDs absent from id, added as founders.} -#' \item{dependencies}{Cycles detected in the pedigree. Must be resolved manually.} +#' \item{dependencies}{Individuals recorded as their own parent, and cycles +#' (individuals that are their own ancestors). Each cycle is listed once, +#' as a path in which every id is a recorded parent of the one before it, +#' back to the first id (e.g. "A -> B -> A"). Must be resolved manually.} #' \item{corrected_pedigree}{Corrected pedigree table.} #' } #' @@ -202,55 +205,96 @@ check_ped <- function(ped.file, } #### check 5: dependencies (cycles) -- reported only #### - detect_all_cycles <- function(data) { - adj_list <- lapply(data$id, function(x) { - row <- data[data$id == x, ] - c(row$male_parent, row$female_parent) - }) - names(adj_list) <- data$id - - dfs <- function(node, visited, rec_stack, path) { - visited[node] <- TRUE - rec_stack[node] <- TRUE - path <- append(path, node) - cycles <- list() - - for (neighbor in adj_list[[node]]) { - if (neighbor %in% names(adj_list)) { - if (!visited[neighbor]) { - cycles <- append(cycles, dfs(neighbor, visited, rec_stack, path)) - } else if (rec_stack[neighbor]) { - cycle_start <- match(neighbor, path) - cycles <- append(cycles, list(path[cycle_start:length(path)])) - } - } + # Individuals that are their own ancestors form a group tied by cycles (a + # strongly connected component of the id -> parent graph). Tarjan's + # algorithm finds each group once, in linear time and without recursion, so + # a cycle is reported once however many individuals it has and however many + # descendants they have. Self-parenting is reported above and skipped here. + cyc_ids <- unique(data$id) + n_ids <- length(cyc_ids) + e_child <- match(c(data$id, data$id), cyc_ids) + e_parent <- match(c(data$male_parent, data$female_parent), cyc_ids) + keep_edge <- !is.na(e_parent) & e_child != e_parent # parents with a row, no self-parenting + e_child <- e_child[keep_edge] + e_parent <- e_parent[keep_edge] + up <- split(e_parent, factor(e_child, levels = seq_len(n_ids))) # parents of each id + + index <- rep(NA_integer_, n_ids) # visit order + low <- integer(n_ids) # lowest visit order reachable + on_stack <- logical(n_ids) + stk_pos <- integer(n_ids) # position of each id in the stack + stk <- integer(n_ids) # ids of the groups being built + path <- integer(n_ids) # current search path (no recursion) + next_up <- integer(n_ids) # next parent to visit, per id + top <- 0L + depth <- 0L + counter <- 0L + groups <- list() + + for (s in seq_len(n_ids)) { + if (!is.na(index[s])) next + todo <- s # id to visit next + while (!is.na(todo) || depth > 0L) { + if (!is.na(todo)) { + # visit: number it and put it on both stacks + counter <- counter + 1L + index[todo] <- counter + low[todo] <- counter + top <- top + 1L + stk[top] <- todo + stk_pos[todo] <- top + on_stack[todo] <- TRUE + depth <- depth + 1L + path[depth] <- todo + next_up[todo] <- 1L + todo <- NA_integer_ } - rec_stack[node] <- FALSE - return(cycles) - } - - visited <- stats::setNames(rep(FALSE, length(adj_list)), names(adj_list)) - rec_stack <- stats::setNames(rep(FALSE, length(adj_list)), names(adj_list)) - all_cycles <- list() - - for (node in names(adj_list)) { - if (!visited[node]) { - node_cycles <- dfs(node, visited, rec_stack, character()) - if (length(node_cycles) > 0) - all_cycles <- append(all_cycles, node_cycles) + v <- path[depth] + if (next_up[v] <= length(up[[v]])) { + w <- up[[v]][next_up[v]] + next_up[v] <- next_up[v] + 1L + if (is.na(index[w])) { + todo <- w + } else if (on_stack[w]) { + low[v] <- min(low[v], index[w]) + } + } else { + # all parents of v done: close v, and its group if v is the root + depth <- depth - 1L + if (depth > 0L) low[path[depth]] <- min(low[path[depth]], low[v]) + if (low[v] == index[v]) { + members <- stk[stk_pos[v]:top] + on_stack[members] <- FALSE + top <- stk_pos[v] - 1L + if (length(members) > 1) groups[[length(groups) + 1]] <- sort(members) + } } } - return(all_cycles) } - - cycles <- detect_all_cycles(data) - if (length(cycles) > 0) { - for (cycle_group in cycles) { - cycle_ids <- unique(unlist(cycle_group)) - errors <- append(errors, - paste("Cycle detected involving IDs:", - paste(cycle_ids, collapse = " -> "))) + + # One line per group, in pedigree order: the loop as a path in which each id + # is a recorded parent of the id before it, back to the first id + groups <- groups[order(vapply(groups, min, integer(1)))] + for (members in groups) { + in_group <- e_child %in% members & e_parent %in% members + g_child <- e_child[in_group] + g_parent <- e_parent[in_group] + loop <- members[1] + v <- g_parent[match(members[1], g_child)] + while (!v %in% loop) { + loop <- c(loop, v) + v <- g_parent[match(v, g_child)] + } + loop <- c(loop[match(v, loop):length(loop)], v) + # a single loop when every id has one parent inside the group + single <- !any(duplicated(g_child[!duplicated(cbind(g_child, g_parent))])) + cycle <- if (single) { + paste(cyc_ids[loop], collapse = " -> ") + } else { + paste0(paste(cyc_ids[members], collapse = ", "), " (more than one loop, e.g. ", + paste(cyc_ids[loop], collapse = " -> "), ")") } + errors <- append(errors, paste("Cycle detected involving IDs:", cycle)) } #### compile findings #### @@ -259,8 +303,8 @@ check_ped <- function(ped.file, conflicting_trios = conflicting_trios, inconsistent_sex_roles = inconsistent_sex_roles, missing_parents = missing_parents, - dependencies = data.frame(dependency = unique(unlist(errors)), - stringsAsFactors = FALSE), + dependencies = data.frame(dependency = as.character(unique(unlist(errors))), + stringsAsFactors = FALSE), # keeps the column when empty corrected_pedigree = data %>% dplyr::select(-row_number) ) diff --git a/docs/BIGpopA_user_guide.Rmd b/docs/BIGpopA_user_guide.Rmd index 291a5ec..ebcae37 100644 --- a/docs/BIGpopA_user_guide.Rmd +++ b/docs/BIGpopA_user_guide.Rmd @@ -25,13 +25,25 @@ knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, fig.align = "center", fig.width = 7, fig.height = 4.5 ) + +# Chunk option fold_output = "label": keep the code visible but collapse its +# printed output behind a clickable label (used for the long verbose reports) +local({ + default_output = knitr::knit_hooks$get("output") + knitr::knit_hooks$set(output = function(x, options) { + out = default_output(x, options) + if (is.null(options$fold_output) || isFALSE(options$fold_output)) return(out) + label = if (is.character(options$fold_output)) options$fold_output else "Show output" + paste0("\n
", label, "\n\n", out, "\n\n
\n") + }) +}) ``` # Introduction BIGpopA is an R package for pedigree quality control, parentage assignment and breed/line composition from SNP genotypes. It works for diploid and polyploid species and is the analytical engine behind the [Familia](https://github.com/Breeding-Insight/Familia) Shiny app. This guide runs the functions in the order a breeding program usually needs them: -1. Convert genotypes to the dosage table used by every function. -> `vcf_to_dosage()`, `ped_to_dosage()` +1. Read genotypes as allele-B dosages; VCF and PLINK files are converted automatically. -> `vcf_to_dosage()`, `ped_to_dosage()` 2. Is the pedigree record clean? -> `check_ped()` 3. Do the recorded parents agree with the genotypes? -> `validate_pedigree()` 4. Who are the parents of an individual with unknown or wrong parents? -> `find_parentage()` @@ -50,10 +62,10 @@ BIGpopA is an R package for pedigree quality control, parentage assignment and b Remarks on the data (violations give misleading results): -- **Genotypes are allele-B dosages** (0 to ploidy), not allele pairs. Convert VCF or PLINK files with the two converters below. +- **Genotypes are allele-B dosages** (0 to the ploidy number), not allele pairs: the number of copies of the counted allele, so 0, 1 or 2 in a diploid and 0 to 4 in a tetraploid. VCF and PLINK files are converted automatically; `vcf_to_dosage()` and `ped_to_dosage()` show the converted table (see [Converting genotype files](#converting-genotype-files)). - **Ploidy is correct.** The Mendelian rules depend on the gametes each parent can produce (ploidy / 2 alleles). A wrong ploidy produces many false failures, so set it to the ploidy of your dosages. Even ploidies use a polysomic gamete-range test (assumes autopolyploid inheritance, so it is conservative for allopolyploids). Odd ploidies (e.g. triploid) fall back to an opposite-homozygote test on homozygous markers only, which has less power. -- **Offspring come from normal sexual reproduction.** Selfed, inbred, doubled-haploid or clonal material can break the assumptions of the trio test (see the Tips section). -- **Marker names and individual ids match exactly** between files (case and spaces matter). +- **Offspring come from normal sexual reproduction.** Selfed, inbred, doubled-haploid or clonal material can break the assumptions of the trio test (see [Tips](#tips)). +- **Individual ids and marker names match exactly** between files (case and spaces matter). Column names are more flexible; see [Input data](#input-data). BIGpopA does not call dosages; use your pipeline first (we recommend [BIGapp](https://github.com/Breeding-Insight/BIGapp)). @@ -68,20 +80,23 @@ library(BIGpopA) # Input data -Every function accepts either a file path or an in-memory object (data.frame / data.table). `founders_file` accepts a file path, a character vector of ids, or a data.frame whose first column holds the ids. This guide builds a small toy dataset that exercises every result category. +BIGpopA works with up to five inputs. Each can be a file path or an in-memory object (data.frame / data.table); genotypes can also be a VCF or a PLINK `.ped` file. This guide builds a small toy dataset that covers the main result categories. -| Object | Format | Used by | -|---|---|---| -| Pedigree | `id`, `male_parent`, `female_parent`. Unknown parent = `0`. | `check_ped()`, `validate_pedigree()` | -| Genotypes | `id` + one column per marker (allele-B dosage, `NA` = missing) | `validate_pedigree()`, `find_parentage()` | -| Parents | `id` and `sex` (`M`, `F`, `A` = either). Without a `sex` column all candidates are treated as `A`. | `find_parentage()` | -| Progeny | `id` column. For *fill_pedigree*: the 3-column pedigree (unknown parent = `0`, `NA` or empty). | `find_parentage()` | -| Founders | File path, character vector or data.frame (first column) | `validate_pedigree()`, `find_parentage()` | +| Input | What it holds | Format | Used by | +|---|---|---|---| +| Pedigree | The recorded parents of each individual, one row per individual | Columns `id`, `male_parent` and `female_parent`. These names are required but matched ignoring case, spaces and dots (`ID` or `Male parent` also work). Unknown parent = `0`. | `check_ped()`, `validate_pedigree()` | +| Genotypes | The marker calls of every genotyped individual, parents and offspring in the same table | `id` + one column per marker with the allele-B dosage (0 to ploidy), missing call = `NA`; or a VCF / PLINK `.ped` file | `validate_pedigree()`, `find_parentage()` | +| Parents | The candidates `find_parentage()` searches | `id` and an optional `sex` (`M`, `F`, or `A` = either). Without `sex` all candidates are `A`. | `find_parentage()` | +| Progeny | The individuals whose parents you want to find | `id`. For `method = "fill_pedigree"`, a pedigree instead (unknown parent = `0`, `NA` or empty). | `find_parentage()` | +| Founders (optional) | The base-generation individuals, whose parents are unknown by design | A vector of ids, a one-column file or a data.frame (ids in the first column) | `validate_pedigree()`, `find_parentage(method = "fill_pedigree")` | + +The same genotype table is used by `validate_pedigree()` and `find_parentage()`. Parents and founders do different jobs. `parents` is the pool of candidates searched by `find_parentage()` and can include any genotyped individual. `founders` only marks the individuals expected to have unknown parents: `validate_pedigree()` reports them as `founders` instead of `missing_both_parents`, and `fill_pedigree` does not search parents for them. In the toy data the four founders are also the four candidate parents. ```{r} library(BIGpopA) -# SNP genotypes coded 0/1/2 (NA = missing); C_low is mostly missing +# Genotypes: allele-B dosage 0/1/2 (diploid), NA = missing call. +# Parents and offspring are in the same table; C_low is mostly missing. genotypes = data.frame( id = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low"), snp01 = c( 0, 0, 2, 2, 0, 2, 1, 0, NA), @@ -99,8 +114,9 @@ genotypes = data.frame( stringsAsFactors = FALSE ) -# pedigree with founders coded 0, a mis-assigned trio (C_bad), -# and a duplicated row (C_missing) +# Pedigree: one row per individual, unknown parents coded 0. It has four +# founders (M1, F1, M2, F2), a mis-assigned trio (C_bad), a sample with mostly +# missing genotypes (C_low) and a duplicated row (C_missing) pedigree = data.frame( id = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low","C_missing","C_missing"), male_parent = c("0", "0", "0", "0", "M1","M2","M1","M2", "M1", "M1", "M1"), @@ -108,7 +124,10 @@ pedigree = data.frame( stringsAsFactors = FALSE ) -# candidate parents (with sex) and progeny to assign +# Founders (optional): individuals whose parents are unknown by design +founders = c("M1", "F1", "M2", "F2") + +# For find_parentage(): candidate parents (with sex) and progeny to assign parents = data.frame( id = c("M1","M2","F1","F2"), sex = c("M", "M", "F", "F"), @@ -119,35 +138,115 @@ progeny = data.frame(id = c("C1","C2","C3","C_bad","C_low","C_missing"), stringsAsFactors = FALSE) pedigree +genotypes ``` -The pedigree has four founders (`M1`, `F1`, `M2`, `F2`, both parents `0`) and several crosses, including a duplicated row (`C_missing`), a mis-assigned trio (`C_bad`) and a sample with mostly missing genotypes (`C_low`). +`C_bad` is recorded as `M2` x `F1` but is really a child of `M1` x `F1`, and `C_missing` is not genotyped. The next sections show how each function handles these cases. ## Converting genotype files -Use these converters when your genotypes come as a VCF or a PLINK file. Both return an `id` column followed by one column per marker, ready for `validate_pedigree()` and `find_parentage()`. +VCF and PLINK files do not need to be converted first: `validate_pedigree()`, `find_parentage()`, `allele_freq_poly()` and `solve_composition_poly()` accept the file path and convert it with their `ploidy`. `vcf_to_dosage()` (VCF, any ploidy) and `ped_to_dosage()` (PLINK `.ped`, diploid) run the same conversion when you want to see or reuse the table. Below are three markers of the toy data for `M1`, `F1` and `C1` in each format. + +```{r, include = FALSE} +# Example files shown in this section (hidden from the guide) +example_dir = tempfile("bigpopa_example_") +dir.create(example_dir) +vcf_file = file.path(example_dir, "example.vcf") +ped_file = file.path(example_dir, "example.ped") +map_file = file.path(example_dir, "example.map") + +writeLines(c( + "##fileformat=VCFv4.3", + '##FORMAT=', + paste("#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO", "FORMAT", + "M1", "F1", "C1", sep = "\t"), + paste("chr1", "100", "snp03", "A", "G", ".", ".", ".", "GT", "0/0", "1/1", "0/1", sep = "\t"), + paste("chr1", "200", "snp05", "C", "T", ".", ".", ".", "GT", "1/1", "0/0", "0/1", sep = "\t"), + paste("chr1", "300", "snp09", "G", "A", ".", ".", ".", "GT", "0/0", "0/1", "0/1", sep = "\t") +), vcf_file) + +writeLines(c("FAM M1 0 0 1 -9 A A T T A A", + "FAM F1 0 0 2 -9 G G C C A G", + "FAM C1 M1 F1 0 -9 A G C T A G"), ped_file) + +writeLines(c("1 snp03 0 100", + "1 snp05 0 200", + "1 snp09 0 300"), map_file) +``` -```{r, eval=FALSE} -# VCF (.vcf or .vcf.gz); ploidy is inferred from the GT field when NULL -geno_vcf = vcf_to_dosage("my_markers.vcf.gz", ploidy = 4) +**example.vcf** (tab-separated; the `GT` calls `0/0`, `0/1` and `1/1` are 0, 1 and 2 copies of the `ALT` allele) -# PLINK .ped (diploid only); a .map file next to it is used for marker names -geno_ped = ped_to_dosage("my_markers.ped") +```{r, echo = FALSE, comment = ""} +cat(readLines(vcf_file), sep = "\n") +``` + +**example.ped** (six pedigree columns, of which only the individual id in column 2 is used, then two alleles per marker) -# format = "matrix" returns individuals x markers, the input of -# allele_freq_poly() and solve_composition_poly() -geno_mat = vcf_to_dosage("my_markers.vcf.gz", ploidy = 4, format = "matrix") +```{r, echo = FALSE, comment = ""} +cat(readLines(ped_file), sep = "\n") ``` -Useful `vcf_to_dosage()` arguments: +**example.map** (one row per marker: chromosome, marker name, genetic position, physical position) -| Argument | Effect | -|---|---| -| `ploidy` | Checked against the ploidy observed in the file; a warning is issued if they disagree | -| `dosage_field` | Read a FORMAT field with dosages (e.g. `"DS"`, `"UD"`) instead of counting alleles in `GT`; `ploidy` must be given | -| `ref_dosage` | Count reference instead of alternate alleles (default `FALSE`) | -| `biallelic_only` | Drop markers with more than one ALT allele (default `TRUE`) | -| `min_marker_call_rate`, `min_sample_call_rate` | Drop markers / individuals with a lower call rate (default 0 = no filtering) | +```{r, echo = FALSE, comment = ""} +cat(readLines(map_file), sep = "\n") +``` + +Converting either file gives the same table: + +```{r, eval = FALSE} +geno_vcf = vcf_to_dosage("example.vcf", verbose = FALSE) # ploidy read from the GT calls +geno_vcf +``` + +```{r, echo = FALSE, eval = requireNamespace("vcfR", quietly = TRUE)} +geno_vcf = vcf_to_dosage(vcf_file, verbose = FALSE) +geno_vcf +``` + +```{r, eval = FALSE} +geno_ped = ped_to_dosage("example.ped", verbose = FALSE) # example.map in the same folder gives the marker names +geno_ped +``` + +```{r, echo = FALSE} +geno_ped = ped_to_dosage(ped_file, verbose = FALSE) +geno_ped +``` + +Both have the layout of `genotypes` and, here, the same values: + +```{r} +genotypes[genotypes$id %in% c("M1", "F1", "C1"), c("id", "snp03", "snp05", "snp09")] +``` + +In a `.ped` file the counted allele is the alphabetically last allele at each marker (`G`, `T` and `G` here). For breed composition, `format = "matrix"` returns the same values with the ids as row names instead of an `id` column: + +```{r, eval = FALSE} +geno_mat = vcf_to_dosage("example.vcf", format = "matrix", verbose = FALSE) +geno_mat +``` + +```{r, echo = FALSE, eval = requireNamespace("vcfR", quietly = TRUE)} +geno_mat = vcf_to_dosage(vcf_file, format = "matrix", verbose = FALSE) +geno_mat +``` + +See `?vcf_to_dosage` and `?ped_to_dosage` for the other options, such as reading a dosage field (`DS`), call-rate filters and multiallelic markers. + +## Missing data + +There are two kinds of missing information, coded differently: + +- **Unknown parent** (pedigree): `0`. In a `fill_pedigree` progeny table `NA` or an empty cell also work. +- **Missing genotype call**: `NA` in a dosage table. The converters turn a VCF call containing `.` (for example `./.`) and a `.ped` allele coded `0` into `NA`. + +Missing calls never count as errors. Each test only uses the markers where every individual it compares has a call (the tests are explained in [How the Mendelian error is calculated](#how-the-mendelian-error-is-calculated)): + +- **Trio test** (`validate_pedigree()`, and `find_parentage()` with `best_pair` or `fill_pedigree`): markers where the offspring and both parents have a call, the *testable* markers. A trio with fewer testable markers than `min_markers` (default 10) is `low_markers`; one with none is `no_data`. +- **Single-parent test** (recommended corrections, and the single-parent methods of `find_parentage()`): markers where the parent and the offspring are both homozygous. If there are none, the parent counts as failing (see `C_low` in [Failed trios and recommended correction](#failed-trios-and-recommended-correction)). +- **Marker summary**: `mismatch_pct` is `NA` for a marker not tested in any trio. +- **Breed composition**: a sample's missing calls are skipped marker by marker. A marker missing in every reference individual of a line has no frequency (`NaN` in the output of `allele_freq_poly()`), which makes `solve_composition_poly()` stop with an error. Find these markers with `rownames(freq)[rowSums(is.na(freq)) > 0]` and remove them before rerunning. # Pedigree cleaning @@ -159,76 +258,134 @@ Useful `vcf_to_dosage()` arguments: | Conflicting trios | Same `id`, different parents in different rows | Disagreeing parents set to `0` and rows collapsed, or flagged (`correct_conflicting_trios = FALSE`) | | Missing parents | A parent is listed but has no row of its own | Always added as a row with unknown parents (`0`) | | Inconsistent sex roles | Same individual used as sire and as dam | Set to `0` as a parent everywhere, or flagged (`correct_inconsistent_sex_roles = FALSE`) | -| Cycles / dependencies | An individual is its own ancestor, or is recorded as its own parent | Flagged only; must be fixed manually, since there is no safe automatic choice | +| Cycles / dependencies | An individual is its own ancestor, or is recorded as its own parent | Flagged only; must be fixed manually, since there is no safe automatic choice. Each cycle is listed once, as a path in which every id is a recorded parent of the one before it (`A -> B -> A`) | + +Unknown parents must be coded `0`. With `verbose = TRUE` the function prints a report of every check: + +```{r} +clean_ped_results = check_ped(pedigree, verbose = TRUE) +``` -Unknown parents must be coded `0`. +The result is a list with one table per issue plus the corrected pedigree. An empty table means the issue was not found: ```{r} -clean_ped_results = check_ped(pedigree, verbose = FALSE) +names(clean_ped_results) +sapply(clean_ped_results, nrow) # number of rows in each table +``` + +| Element | In the toy data | +|---|---| +| `exact_duplicates` | Both copies of the repeated `C_missing` row; `row_number` is their position in the input | +| `conflicting_trios` | Empty: no id has two different sets of parents | +| `inconsistent_sex_roles` | Empty: no individual is used as both male and female parent | +| `missing_parents` | Empty: every parent has its own row | +| `dependencies` | Empty: no cycles | +| `corrected_pedigree` | The pedigree without the duplicate, ready for `validate_pedigree()` | + +```{r} +clean_ped_results$exact_duplicates # the corrected, analysis-ready pedigree clean_ped = clean_ped_results$corrected_pedigree clean_ped ``` -The result also holds one table per issue, for example `clean_ped_results$exact_duplicates` (here the repeated `C_missing` row) and `clean_ped_results$missing_parents`. An empty table means the issue was not found. Set `verbose = TRUE` to print a full report to the console. - # Pedigree validation -`validate_pedigree()` compares every trio (offspring, male parent, female parent) with the genotypes and computes a Mendelian error rate. A marker is *testable* when the offspring and both parents have a call. +`validate_pedigree()` tests every trio (offspring, male parent, female parent) of the pedigree against the genotypes. It returns the status of each trio, a recommended correction when a trio fails, a corrected pedigree and a per-marker summary. Here the founders are passed with `founders_file`, and `marker_trio_table = TRUE` adds the marker-by-trio table used below. Click the label under the code to see the console report. -```{r} +```{r, fold_output = "Show the report printed with verbose = TRUE"} ped_validate_results = validate_pedigree(clean_ped, genotypes, - verbose = FALSE, plot_results = TRUE) + founders_file = founders, + marker_trio_table = TRUE, + verbose = TRUE, plot_results = TRUE) +``` + +The result is a list. `full_results` has one row per trio, and `pass`, `fail`, `low_markers`, `no_genotype_data`, `founders` and `missing_parents` hold the same rows split by status. `corrected_pedigree`, `marker_summary` and `marker_trio_table` are described below, and `plot` is the histogram above. + +```{r} +names(ped_validate_results) +sapply(Filter(is.data.frame, ped_validate_results), nrow) # rows in each table ped_validate_report = ped_validate_results$full_results ped_validate_report ``` -## Parameters +In `full_results`, `trio_mendelian_error_pct` is the trio error over `trio_markers_tested` testable markers, and `male_parent_hom_error_pct` and `female_parent_hom_error_pct` are the single-parent errors, computed only for failed and low-marker trios. The thresholds used here are the defaults; they and the other settings are arguments of the function (see `?validate_pedigree`). + +## How the Mendelian error is calculated + +The trio error uses **every marker where the offspring and both parents have a call**, heterozygous or homozygous, not only the homozygous ones. At each testable marker the question is whether the two parents can produce the offspring. Each parent passes on a gamete with half of its alleles (ploidy / 2), so a parent with dosage *g* passes between max(0, *g* − ploidy / 2) and min(*g*, ploidy / 2) copies of allele B. The offspring dosage must lie between the sum of the two parents' minimums and the sum of their maximums; otherwise the marker is a Mendelian error. For a diploid this gives the familiar rules: -| Argument | Default | Effect | +| Male parent | Female parent | Possible offspring | |---|---|---| -| `ploidy` | 2 | Ploidy of the dosages | -| `trio_error_threshold` | 5 | Maximum % of testable markers with a Mendelian error for a trio to *pass* | -| `min_markers` | 10 | Minimum testable markers needed to judge a trio; fewer gives `low_markers` | -| `single_parent_error_threshold` | 2 | Maximum homozygous mismatch % for a single parent to pass the single-parent check | -| `founders_file` | `NULL` | File path, vector of ids or data.frame; ids listed here with both parents `0` are reported as `founders` | -| `marker_trio_table` | `FALSE` | Also return the marker x trio table (see below) | -| `plot_results` | `TRUE` | Draw a histogram of trio error percentages | +| 0 | 0 | 0 | +| 0 | 1 | 0 or 1 | +| 0 | 2 | 1 | +| 1 | 1 | 0, 1 or 2 | +| 1 | 2 | 1 or 2 | +| 2 | 2 | 2 | + +For a tetraploid, a parent with dosage 3 passes 1 or 2 copies, so a 3 x 0 cross can only give 1 or 2. The trio error is then + +> trio error (%) = 100 × markers with an error / testable markers + +and the trio passes when it is at or below `trio_error_threshold` (default 5%). For odd ploidy (e.g. triploid) gametes are unbalanced, so the trio test falls back to the opposite-homozygote check of the single-parent test (see [Failed trios and recommended correction](#failed-trios-and-recommended-correction)): a marker is an error when the offspring is homozygous and a homozygous parent carries the other allele. + +Here is the trio of `C_bad` (recorded as `M2` x `F1`) marker by marker, with the verdict from `marker_trio_table`: + +```{r} +c_bad_trio = data.frame( + marker = ped_validate_results$marker_trio_table$marker, + M2 = unlist(genotypes[genotypes$id == "M2", -1]), + F1 = unlist(genotypes[genotypes$id == "F1", -1]), + C_bad = unlist(genotypes[genotypes$id == "C_bad", -1]), + verdict = ped_validate_results$marker_trio_table$C_bad_M2_F1, + row.names = NULL +) +c_bad_trio +``` + +`r sum(c_bad_trio$verdict != "match")` of the `r ped_validate_report$trio_markers_tested[ped_validate_report$id == "C_bad"]` testable markers are errors, a trio error of `r ped_validate_report$trio_mendelian_error_pct[ped_validate_report$id == "C_bad"]`%, so the trio fails. At `snp01`, for example, `M2` is 2 and `F1` is 0, so every child must be 1, but `C_bad` is 0. Heterozygous markers count too: at `snp10` the cross 2 x 1 can only give 1 or 2, and `C_bad` is 0. All errors are labelled `mismatch_male`, meaning `M2` cannot be the parent whatever `F1` contributes; the single-parent test below confirms this. ## Status of each trio | Status | Meaning | |---|---| -| `pass` | Mendelian error at or below `trio_error_threshold` | -| `fail` | Error above the threshold | +| `pass` | Trio error at or below `trio_error_threshold` | +| `fail` | Trio error above the threshold | | `low_markers` | Fewer than `min_markers` testable markers | -| `no_genotype_data` | At least one of the three individuals is not in the genotype file | -| `founders` | Both parents `0` and the id is in `founders_file` | -| `missing_both_parents`, `missing_male_parent`, `missing_female_parent` | One or both parents are `0` | -| `no_data` | No testable marker at all | +| `no_genotype_data` | At least one of the three individuals is not in the genotype table | +| `founders` | Both parents `0` and the id is listed in `founders_file` | +| `missing_both_parents`, `missing_male_parent`, `missing_female_parent` | One or both parents are `0`, and the individual is not a listed founder | +| `no_data` | No testable marker at all (only in `full_results`) | -In the toy data `C1`, `C2` and `C3` pass, `C_bad` fails, `C_low` is `low_markers` and `C_missing` has `no_genotype_data`. +In the toy data `C1`, `C2` and `C3` pass, `C_bad` fails, `C_low` is `low_markers` and `C_missing` has `no_genotype_data`. `M1`, `F1`, `M2` and `F2` are `founders` because they are listed in `founders`; without `founders_file` they would be `missing_both_parents`. ## Failed trios and recommended correction -For a failed trio the function checks each parent alone, using only markers where the parent and the offspring are both homozygous (an offspring homozygous for one allele cannot have a parent homozygous for the other). The `recommended_correction` column gives the result: +A failed trio shows that the pair is wrong, but not which parent. Each parent is therefore tested alone with the **single-parent test**, which uses only the markers where the parent and the offspring are both homozygous. A parent homozygous for one allele passes that allele to every child, so a child homozygous for the other allele (0 vs 2 in a diploid) cannot be its offspring. Heterozygous markers are left out because, without the other parent, they rarely rule a parent out. + +> single-parent error (%) = 100 × markers where parent and offspring are homozygous for different alleles / markers where both are homozygous + +A parent fails when this error is above `single_parent_error_threshold` (default 2%), or when there is no marker to compare. For `C_bad`, `M2` has a single-parent error of `r ped_validate_report$male_parent_hom_error_pct[ped_validate_report$id == "C_bad"]`% and `F1` of `r ped_validate_report$female_parent_hom_error_pct[ped_validate_report$id == "C_bad"]`%, so the recommendation is `remove_male_parent`. The `recommended_correction` column gives the result for every trio: | `recommended_correction` | Interpretation | |---|---| | `none` | Trio passes, or it is a founder, has a missing parent, has no genotype data or has no testable marker | -| `remove_male_parent` | The trio fails and only the male parent fails the single-parent check | -| `remove_female_parent` | The trio fails and only the female parent fails the single-parent check | -| `remove_both` | The trio fails and both parents fail the single-parent check | -| `unresolved` | The trio fails but neither parent fails the single-parent check, for example an offspring recorded as its own parent; the pedigree is probably wrong but the faulty parent cannot be identified | +| `remove_male_parent` | The trio fails and only the male parent fails the single-parent test | +| `remove_female_parent` | The trio fails and only the female parent fails the single-parent test | +| `remove_both` | The trio fails and both parents fail the single-parent test | +| `unresolved` | The trio fails but neither parent fails the single-parent test, for example an offspring recorded as its own parent; the pedigree is probably wrong but the faulty parent cannot be identified | | `low_markers_*` | The trio has fewer than `min_markers` testable markers; the suffix is the decision the same rules would give (`low_markers_keep_both` means the error is within the threshold and both parents pass) | +`C_low` shows the weak point of having few markers. Its trio error is `r ped_validate_report$trio_mendelian_error_pct[ped_validate_report$id == "C_low"]`%, but `F1` is heterozygous at all `r ped_validate_report$trio_markers_tested[ped_validate_report$id == "C_low"]` of its testable markers, so the single-parent test has nothing to compare for `F1`. `F1` then counts as failing, and the recommendation is `low_markers_remove_female_parent`. Treat `low_markers_*` recommendations with caution, and genotype these samples again if possible. + ```{r} ped_validate_results$corrected_pedigree ``` -`corrected_pedigree` sets the parent(s) named by the recommendation to `0`. Run `find_parentage()` with `method = "fill_pedigree"` afterwards to look for the right parent. +`corrected_pedigree` sets the parent(s) named by the recommendation to `0`, here the male parent of `C_bad` and the female parent of `C_low`. `unresolved` trios are left unchanged for manual review. Run `find_parentage()` with `method = "fill_pedigree"` on the corrected pedigree to look for the right parents (see [Filling in missing parents of a pedigree](#filling-in-missing-parents-of-a-pedigree)). ## Marker summary @@ -240,12 +397,10 @@ head(ped_validate_results$marker_summary) ## Marker x trio table -Set `marker_trio_table = TRUE` to see, for each marker and each trio, where an error comes from. Columns are named `id_male_female`. +With `marker_trio_table = TRUE` (as in the run above), the result includes the verdict of every marker for every trio with both parents known and genotyped. Columns are named `id_male_female`; the `C_bad_M2_F1` column was used above. ```{r} -mt_results = validate_pedigree(clean_ped, genotypes, verbose = FALSE, - plot_results = FALSE, marker_trio_table = TRUE) -mt_results$marker_trio_table +ped_validate_results$marker_trio_table ``` | Label | Meaning | @@ -257,74 +412,71 @@ mt_results$marker_trio_table | `missing_male`, `missing_female`, `missing_progeny` | That individual has no call; every missing individual is named (for example `missing_progeny_male_female`) | | `untestable` | No call is missing, but the marker cannot be tested (odd ploidy, offspring not homozygous) | -# Missing data +# Parentage assignment -Missing marker calls (`NA`) are handled marker by marker and never count as errors: +When parentage is unknown, `find_parentage()` tests every candidate in `parents` against each individual in `progeny` and reports the best one. How candidates are scored depends on the method; the score is reported as `mendelian_error_pct`, and lower is better: -- **validate_pedigree().** A marker is *testable* for a trio only if the offspring and both parents have a call (odd ploidy: the offspring is homozygous and at least one parent has a call and is homozygous). The error % is mismatches / testable markers. Trios with fewer testable markers than `min_markers` are labeled `low_markers`; trios with none are `no_data` (only in `full_results`). -- **Single-parent check.** Uses only markers where both parent and offspring are homozygous. If there are none, the parent counts as failing. -- **Marker summary.** `mismatch_pct` is `NA` for markers not tested in any trio. -- **find_parentage().** Same testable-marker logic per candidate; a candidate with no testable markers gets `markers_tested = 0` and `low_markers`. -- **allele_freq_poly().** A marker missing in every reference individual of a line has no frequency; `solve_composition_poly()` stops with an error listing these markers. Remove them and rerun. +| Method | Candidates tested | Score | +|---|---|---| +| `best_pair` (default) | Every male x female pair | Trio error, as in `validate_pedigree()`: all markers where the three individuals have a call | +| `best_male_parent`, `best_female_parent` | Each candidate of that sex alone | Single-parent error: only markers homozygous in both candidate and progeny | +| `best_match` | Each candidate of either sex alone | Single-parent error, as above | +| `fill_pedigree` | The missing parent(s) of a pedigree; a known parent is kept fixed | Trio error (see [Filling in missing parents of a pedigree](#filling-in-missing-parents-of-a-pedigree)) | -# Parentage assignment +Only the single-parent methods rely on homozygous markers. Prefer `best_pair` when both parents are among the candidates: it uses all markers and both parents together, so it separates candidates much better. Use the single-parent methods when the other parent is not genotyped. -When parentage is unknown, `find_parentage()` searches a pool of candidate parents for the best assignment to each progeny. +The best candidate is always reported, with a status: `pass` (score at or below `error_threshold`, default 5%), `high_error` (above it: no candidate fits well, so the true parent may be missing from `parents`), `low_markers` (fewer than `min_markers` markers tested) or, for `fill_pedigree`, `known_parent_no_genotype`. See `?find_parentage` for all options. -```{r} +```{r, fold_output = "Show the report printed with verbose = TRUE"} find_parentage_results = find_parentage(genotypes, parents, progeny, method = "best_pair", - verbose = FALSE, plot_results = TRUE) + verbose = TRUE, plot_results = TRUE) +``` + +The result is a list like the one from `validate_pedigree()`: `full_results` has one row per progeny, `pass`, `high_error` and `low_markers` hold the same rows split by status, and `plot` is the histogram above. + +```{r} +names(find_parentage_results) +sapply(Filter(is.data.frame, find_parentage_results), nrow) # rows in each table find_parentage_report = find_parentage_results$full_results find_parentage_report ``` -## Parameters - -| Argument | Default | Effect | -|---|---|---| -| `method` | `"best_pair"` | The search (see below) | -| `ploidy` | 2 | Ploidy of the dosages | -| `error_threshold` | 5 | Error % above which the best candidate gets `high_error` | -| `min_markers` | 10 | If the best candidate has fewer testable markers, the status is `low_markers` | -| `top_n` | 1 | Report the best `top_n` candidates per progeny | -| `show_ties` | `TRUE` | Report `tied_candidates` | -| `allow_parent_selfing` | `FALSE` | Allow the same individual as male and female parent | -| `exclude_self_match` | `TRUE` | Do not let a progeny match itself as parent | -| `founders_file` | `NULL` | File path, vector of ids or data.frame; founders are skipped by `fill_pedigree` | - -| Method | Searches | -|---|---| -| `best_pair` | Male x female pair with the lowest Mendelian error (needs both sexes in the parents table) | -| `best_male_parent`, `best_female_parent` | The single candidate of that sex with the lowest homozygous mismatch (compares only markers homozygous in both candidate and offspring); use when the other parent is known or unavailable | -| `best_match` | Same homozygous mismatch test over candidates of either sex | -| `fill_pedigree` | Missing parents of a 3-column pedigree (see below) | - -The best candidate is always reported. The status is `pass` (error at or below `error_threshold`), `high_error` (above it, so treat the assignment with caution), `low_markers`, or `known_parent_no_genotype`. In the toy data `C_bad`, which failed validation, is recovered to its true pair (`M1`/`F1`) at 0% error, while `C_low` remains `low_markers`. +Each row gives the best pair (`male_parent`, `female_parent`), its error (`mendelian_error_pct`), the number of markers tested, the status and `tied_candidates`, the number of pairs that share the best error. `C1`, `C2` and `C3` get their recorded parents at 0% error, and `C_bad`, which failed validation, is recovered to its true pair (`M1` x `F1`). `C_low` also gets `M1` x `F1`, but on only `r find_parentage_report$markers_tested[find_parentage_report$id == "C_low"]` markers, so it stays `low_markers`. `C_missing` is not in the genotype table, so it is dropped with a warning. ## Top N and ties -`tied_candidates` is the number of candidates (pairs for `best_pair`) that share the lowest error of a progeny; values above 1 mean the best assignment is ambiguous. This is common with few markers or with clones. With `top_n > 1`, rank 1 stays in the base columns and ranks 2 to `top_n` are added as suffix columns (`male_parent_2`, `female_parent_2`, `mendelian_error_pct_2`). +`tied_candidates` above 1 means several candidates share the best error, so the assignment is ambiguous; this is common with few markers or with clones. To see how close the alternatives are, set `top_n`: ```{r} -find_parentage(genotypes, parents, progeny, method = "best_pair", top_n = 3, - verbose = FALSE, plot_results = FALSE)$full_results +top3_results = find_parentage(genotypes, parents, progeny, method = "best_pair", + top_n = 3, verbose = FALSE, plot_results = FALSE) +names(top3_results) + +top3 = top3_results$full_results +top3 ``` +The list has the same elements; only `full_results` gains columns. Rank 1 stays in the base columns (`male_parent`, `female_parent`, `mendelian_error_pct`, `markers_tested`), and ranks 2 and 3 are added with the suffixes `_2` and `_3` (`male_parent_2`, `female_parent_2`, `mendelian_error_pct_2`, `markers_tested_2`, and so on). `status` and `tied_candidates` refer to rank 1. A large gap between ranks means a clear assignment: for `C_bad` the best pair has `r top3$mendelian_error_pct[top3$id == "C_bad"]`% error and the second best (`r top3$male_parent_2[top3$id == "C_bad"]` x `r top3$female_parent_2[top3$id == "C_bad"]`) `r top3$mendelian_error_pct_2[top3$id == "C_bad"]`%. With the default `top_n = 1`, tied best pairs are added in the same suffix columns (`show_ties = TRUE`). + ## Filling in missing parents of a pedigree -With `method = "fill_pedigree"`, `progeny` is a pedigree with `id`, `male_parent` and `female_parent` (for example `validate_pedigree()$corrected_pedigree`). Parents coded `0`, `NA` or empty are unknown. A known male parent searches only females, a known female parent only males, and rows with both unknown get the full pair search; rows with both parents known are skipped. The parents table (with sex `M`, `F` or `A`) is still required. The output adds `search_mode` and `known_parent_error_pct`. +With `method = "fill_pedigree"`, `progeny` is a pedigree (`id`, `male_parent`, `female_parent`), typically the `corrected_pedigree` from `validate_pedigree()`. Parents coded `0`, `NA` or empty are unknown. Each row is handled on its own: a known male parent searches only females, a known female parent only males, and a row with both parents unknown gets the full pair search, unless it is a listed founder. Rows with both parents known are skipped. The `parents` table (with sex) is still required. ```{r} -ped_to_fill = data.frame(id = c("C1", "C2"), - male_parent = c("M1", "0"), - female_parent = c("0", "0")) - -find_parentage(genotypes, parents, ped_to_fill, method = "fill_pedigree", - verbose = FALSE, plot_results = FALSE)$full_results +filled_results = find_parentage(genotypes, parents, + ped_validate_results$corrected_pedigree, + method = "fill_pedigree", + founders_file = founders, + verbose = FALSE, plot_results = FALSE) +names(filled_results) + +filled_results$full_results ``` +The list has the same elements as before. Only the rows with a parent removed by validation are searched: the founders are skipped because they are in `founders`, and the complete trios because both parents are known. `full_results` gains two columns: `search_mode` (`known_male`, `known_female`, or `pair` when both parents are unknown) and `known_parent_error_pct`, the single-parent error of the parent that was kept (a high value suggests that it is wrong too). `C_bad` kept `F1`, so only males were tested, and `M1` is found at 0%. `C_low` kept `M1` and gets `F1` back, but on too few markers, so it stays `low_markers`. + # Breed and line composition BIGpopA estimates genome-wide breed or line composition for diploid and polyploid populations with a breedTools-style quadratic programming approach (Funkhouser et al. 2017). The workflow has three steps: compute reference-line allele frequencies, estimate the composition of each sample, and assign each sample to a line. This example mirrors an autotetraploid analysis (`ploidy = 4`) that assigns samples to two cultivar reference panels, Jewel and Draper. @@ -371,7 +523,7 @@ prediction = as.data.frame(solve_composition_poly(validation, freq, ploidy = 4)) prediction ``` -The model returns an R2 value, but it is expected to be low because a discrete allele dosage is regressed on a continuous allele frequency. +`R2` measures how well the estimated mixture of reference lines explains a sample's genotypes. The sample's dosages are divided by the ploidy (0, 0.25, ..., 1 for a tetraploid); the predicted value at each marker is the allele frequency of the reference lines weighted by the estimated proportions; and `R2` is the squared correlation between the two across the sample's non-missing markers. It describes the fit, not the accuracy of the proportions. It is usually well below 1 even for a pure individual, because each genotype is one of a few discrete values while the prediction is a frequency, so every marker at which the line is polymorphic adds unexplained variation. Compare `R2` between samples analysed with the same reference panel and markers: a sample with a much lower `R2` than the rest may come from a population that is not in the reference panel, or have poor-quality genotypes. `assign_composition_line()` turns the result into the table used by Familia: the estimated percentage of each panel and a `Predicted line`. A sample is assigned to its highest-proportion line only if that proportion is at least `assign_threshold` (default 50); otherwise it is labelled `Undetermined`. With `assign_threshold = 0` every sample is assigned. @@ -385,18 +537,23 @@ pred_results # Complete workflow ```{r, eval=FALSE} -# genotypes +# genotypes: a dosage table, or pass the VCF / .ped path straight to each function genotypes = vcf_to_dosage("my_markers.vcf.gz", ploidy = 4) # pedigree quality control clean_ped = check_ped(pedigree)$corrected_pedigree -validated = validate_pedigree(clean_ped, genotypes, ploidy = 4) -assigned = find_parentage(genotypes, parents, progeny, method = "fill_pedigree", - ploidy = 4) +validated = validate_pedigree(clean_ped, genotypes, founders_file = founders, ploidy = 4) + +# look for the parents removed by validation (and any other missing parent) +filled = find_parentage(genotypes, parents, validated$corrected_pedigree, + method = "fill_pedigree", founders_file = founders, ploidy = 4) + +# or assign parents to new progeny +assigned = find_parentage(genotypes, parents, progeny, method = "best_pair", ploidy = 4) # breed/line composition -freq = allele_freq_poly(reference, ref_ids, ploidy = 4) -prediction = solve_composition_poly(validation, freq, ploidy = 4) +freq = allele_freq_poly(reference, ref_ids, ploidy = 4) +prediction = solve_composition_poly(validation, freq, ploidy = 4) composition_table = assign_composition_line(prediction, assign_threshold = 50) ``` @@ -404,6 +561,7 @@ composition_table = assign_composition_line(prediction, assign_threshold = 50) - Typical workflow: `check_ped()` -> `validate_pedigree()` -> `find_parentage()` (`fill_pedigree`) on the corrected pedigree. - Make sure you are using the right ploidy. A wrong ploidy gives many false failures. +- List your founders in `founders_file`, so their unknown parents are not reported as missing. - Selfed or inbred offspring can fail the trio test; check the single-parent errors and the `unresolved` flag. - Ties are common with few markers or with clones; always check `tied_candidates`. - Check `marker_summary` before blaming individuals; remove bad markers and rerun. @@ -412,7 +570,7 @@ composition_table = assign_composition_line(prediction, assign_threshold = 50) ## BIGpopA -Chinchilla-Vargas J, Sandercock AM (2026). *BIGpopA: Pedigree Validation Genetic Composition of Diploids & Polyploids*. R package version 2.2.0. CRAN: +Chinchilla-Vargas J, Sandercock AM (2026). *BIGpopA: Pedigree Validation Genetic Composition of Diploids & Polyploids*. R package version `r as.character(packageVersion("BIGpopA"))`. CRAN: Run `citation("BIGpopA")` in R for the BibTeX entry. diff --git a/docs/BIGpopA_user_guide.html b/docs/BIGpopA_user_guide.html index d45ce52..dcc6840 100644 --- a/docs/BIGpopA_user_guide.html +++ b/docs/BIGpopA_user_guide.html @@ -11,7 +11,7 @@ - + BIGpopA Tutorial @@ -3443,7 +3443,7 @@

BIGpopA Tutorial

Pedigree cleaning, validation, parentage assignment and breed composition from SNP genotype data

Josue Chinchilla-Vargas and Alexander Sandercock

-

2026-10-12

+

2026-10-26

@@ -3456,8 +3456,9 @@

1 Introduction

app. This guide runs the functions in the order a breeding program usually needs them:

    -
  1. Convert genotypes to the dosage table used by every function. -> -vcf_to_dosage(), ped_to_dosage()
    2. +
  2. Read genotypes as allele-B dosages; VCF and PLINK files are +converted automatically. -> vcf_to_dosage(), +ped_to_dosage()
  3. Is the pedigree record clean? -> check_ped()
  4. Do the recorded parents agree with the genotypes? -> validate_pedigree()
    5. @@ -3525,9 +3526,11 @@

    1 Introduction

    Remarks on the data (violations give misleading results):

      -
    • Genotypes are allele-B dosages (0 to ploidy), not -allele pairs. Convert VCF or PLINK files with the two converters -below.
      • +
    • Genotypes are allele-B dosages (0 to the ploidy +number), not allele pairs: the number of copies of the counted allele, +so 0, 1 or 2 in a diploid and 0 to 4 in a tetraploid. VCF and PLINK +files are converted automatically; vcf_to_dosage() and +ped_to_dosage() show the converted table (see Converting genotype files).
    • Ploidy is correct. The Mendelian rules depend on the gametes each parent can produce (ploidy / 2 alleles). A wrong ploidy produces many false failures, so set it to the ploidy of your dosages. @@ -3537,9 +3540,10 @@

      1 Introduction

      markers only, which has less power.
    • Offspring come from normal sexual reproduction. Selfed, inbred, doubled-haploid or clonal material can break the -assumptions of the trio test (see the Tips section).
      • -
    • Marker names and individual ids match exactly -between files (case and spaces matter).
      • +assumptions of the trio test (see Tips). +
    • Individual ids and marker names match exactly +between files (case and spaces matter). Column names are more flexible; +see Input data.

    BIGpopA does not call dosages; use your pipeline first (we recommend BIGapp).

    @@ -3553,20 +3557,21 @@

    2 Install

    3 Input data

    -

    Every function accepts either a file path or an in-memory object -(data.frame / data.table). founders_file accepts a file -path, a character vector of ids, or a data.frame whose first column -holds the ids. This guide builds a small toy dataset that exercises -every result category.

    +

    BIGpopA works with up to five inputs. Each can be a file path or an +in-memory object (data.frame / data.table); genotypes can also be a VCF +or a PLINK .ped file. This guide builds a small toy dataset +that covers the main result categories.

    ---++++ - + + @@ -3574,145 +3579,228 @@

    3 Input data

    - + + - + + - + + - + - - - + + + +
    ObjectInputWhat it holds Format Used by
    Pedigreeid, male_parent, -female_parent. Unknown parent = 0.The recorded parents of each individual, one row per individualColumns id, male_parent and +female_parent. These names are required but matched +ignoring case, spaces and dots (ID or +Male parent also work). Unknown parent = +0. check_ped(), validate_pedigree()
    Genotypesid + one column per marker (allele-B dosage, -NA = missing)The marker calls of every genotyped individual, parents and +offspring in the same tableid + one column per marker with the allele-B dosage (0 +to ploidy), missing call = NA; or a VCF / PLINK +.ped file validate_pedigree(), find_parentage()
    Parentsid and sex (M, -F, A = either). Without a sex -column all candidates are treated as A.The candidates find_parentage() searchesid and an optional sex (M, +F, or A = either). Without sex +all candidates are A. find_parentage()
    Progenyid column. For fill_pedigree: the 3-column -pedigree (unknown parent = 0, NA or +The individuals whose parents you want to findid. For method = "fill_pedigree", a +pedigree instead (unknown parent = 0, NA or empty). find_parentage()
    FoundersFile path, character vector or data.frame (first column)validate_pedigree(), find_parentage()Founders (optional)The base-generation individuals, whose parents are unknown by +designA vector of ids, a one-column file or a data.frame (ids in the first +column)validate_pedigree(), +find_parentage(method = "fill_pedigree")
    +

    The same genotype table is used by validate_pedigree() +and find_parentage(). Parents and founders do different +jobs. parents is the pool of candidates searched by +find_parentage() and can include any genotyped individual. +founders only marks the individuals expected to have +unknown parents: validate_pedigree() reports them as +founders instead of missing_both_parents, and +fill_pedigree does not search parents for them. In the toy +data the four founders are also the four candidate parents.

    library(BIGpopA)
 
-# SNP genotypes coded 0/1/2 (NA = missing); C_low is mostly missing
-genotypes = data.frame(
-  id    = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low"),
-  snp01 = c(  0,   0,   2,   2,   0,   2,   1,     0,     NA),
-  snp02 = c(  0,   0,   2,   0,   0,   1,   0,     0,     NA),
-  snp03 = c(  0,   2,   0,   2,   1,   1,   1,     1,     NA),
-  snp04 = c(  0,   2,   0,   0,   1,   0,   0,     1,     NA),
-  snp05 = c(  2,   0,   1,   1,   1,   1,   1,     1,     NA),
-  snp06 = c(  2,   0,   1,   0,   1,   1,   1,     1,     NA),
-  snp07 = c(  2,   2,   0,   1,   2,   0,   2,     2,     NA),
-  snp08 = c(  2,   2,   2,   2,   2,   2,   2,     2,     NA),
-  snp09 = c(  0,   1,   0,   2,   1,   1,   1,     1,      0),
-  snp10 = c(  0,   1,   2,   0,   0,   1,   0,     0,      0),
-  snp11 = c(  2,   1,   1,   2,   1,   2,   2,     1,      2),
-  snp12 = c(  2,   1,   0,   1,   2,   1,   1,     2,      2),
-  stringsAsFactors = FALSE
-)
-
-# pedigree with founders coded 0, a mis-assigned trio (C_bad),
-# and a duplicated row (C_missing)
-pedigree = data.frame(
-  id            = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low","C_missing","C_missing"),
-  male_parent   = c("0", "0", "0", "0", "M1","M2","M1","M2",   "M1",   "M1",       "M1"),
-  female_parent = c("0", "0", "0", "0", "F1","F2","F2","F1",   "F1",   "F1",       "F1"),
-  stringsAsFactors = FALSE
-)
-
-# candidate parents (with sex) and progeny to assign
-parents = data.frame(
-  id  = c("M1","M2","F1","F2"),
-  sex = c("M", "M", "F", "F"),
-  stringsAsFactors = FALSE
-)
-
-progeny = data.frame(id = c("C1","C2","C3","C_bad","C_low","C_missing"),
-                     stringsAsFactors = FALSE)
-
-pedigree
    +# Genotypes: allele-B dosage 0/1/2 (diploid), NA = missing call. +# Parents and offspring are in the same table; C_low is mostly missing. +genotypes = data.frame( + id = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low"), + snp01 = c( 0, 0, 2, 2, 0, 2, 1, 0, NA), + snp02 = c( 0, 0, 2, 0, 0, 1, 0, 0, NA), + snp03 = c( 0, 2, 0, 2, 1, 1, 1, 1, NA), + snp04 = c( 0, 2, 0, 0, 1, 0, 0, 1, NA), + snp05 = c( 2, 0, 1, 1, 1, 1, 1, 1, NA), + snp06 = c( 2, 0, 1, 0, 1, 1, 1, 1, NA), + snp07 = c( 2, 2, 0, 1, 2, 0, 2, 2, NA), + snp08 = c( 2, 2, 2, 2, 2, 2, 2, 2, NA), + snp09 = c( 0, 1, 0, 2, 1, 1, 1, 1, 0), + snp10 = c( 0, 1, 2, 0, 0, 1, 0, 0, 0), + snp11 = c( 2, 1, 1, 2, 1, 2, 2, 1, 2), + snp12 = c( 2, 1, 0, 1, 2, 1, 1, 2, 2), + stringsAsFactors = FALSE +) + +# Pedigree: one row per individual, unknown parents coded 0. It has four +# founders (M1, F1, M2, F2), a mis-assigned trio (C_bad), a sample with mostly +# missing genotypes (C_low) and a duplicated row (C_missing) +pedigree = data.frame( + id = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low","C_missing","C_missing"), + male_parent = c("0", "0", "0", "0", "M1","M2","M1","M2", "M1", "M1", "M1"), + female_parent = c("0", "0", "0", "0", "F1","F2","F2","F1", "F1", "F1", "F1"), + stringsAsFactors = FALSE +) + +# Founders (optional): individuals whose parents are unknown by design +founders = c("M1", "F1", "M2", "F2") + +# For find_parentage(): candidate parents (with sex) and progeny to assign +parents = data.frame( + id = c("M1","M2","F1","F2"), + sex = c("M", "M", "F", "F"), + stringsAsFactors = FALSE +) + +progeny = data.frame(id = c("C1","C2","C3","C_bad","C_low","C_missing"), + stringsAsFactors = FALSE) + +pedigree
    -

    The pedigree has four founders (M1, F1, -M2, F2, both parents 0) and -several crosses, including a duplicated row (C_missing), a -mis-assigned trio (C_bad) and a sample with mostly missing -genotypes (C_low).

    +
    genotypes
    +
    + +
    +

    C_bad is recorded as M2 x F1 +but is really a child of M1 x F1, and +C_missing is not genotyped. The next sections show how each +function handles these cases.

    3.1 Converting genotype files

    -

    Use these converters when your genotypes come as a VCF or a PLINK -file. Both return an id column followed by one column per -marker, ready for validate_pedigree() and -find_parentage().

    -
    # VCF (.vcf or .vcf.gz); ploidy is inferred from the GT field when NULL
-geno_vcf = vcf_to_dosage("my_markers.vcf.gz", ploidy = 4)
-
-# PLINK .ped (diploid only); a .map file next to it is used for marker names
-geno_ped = ped_to_dosage("my_markers.ped")
-
-# format = "matrix" returns individuals x markers, the input of
-# allele_freq_poly() and solve_composition_poly()
-geno_mat = vcf_to_dosage("my_markers.vcf.gz", ploidy = 4, format = "matrix")
    -

    Useful vcf_to_dosage() arguments:

    - ---- - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ArgumentEffect
    ploidyChecked against the ploidy observed in the file; a warning is issued -if they disagree
    dosage_fieldRead a FORMAT field with dosages (e.g. "DS", -"UD") instead of counting alleles in GT; -ploidy must be given
    ref_dosageCount reference instead of alternate alleles (default -FALSE)
    biallelic_onlyDrop markers with more than one ALT allele (default -TRUE)
    min_marker_call_rate, -min_sample_call_rateDrop markers / individuals with a lower call rate (default 0 = no -filtering)
    +

    VCF and PLINK files do not need to be converted first: +validate_pedigree(), find_parentage(), +allele_freq_poly() and +solve_composition_poly() accept the file path and convert +it with their ploidy. vcf_to_dosage() (VCF, +any ploidy) and ped_to_dosage() (PLINK .ped, +diploid) run the same conversion when you want to see or reuse the +table. Below are three markers of the toy data for M1, +F1 and C1 in each format.

    +

    example.vcf (tab-separated; the GT +calls 0/0, 0/1 and 1/1 are 0, 1 +and 2 copies of the ALT allele)

    +
    ##fileformat=VCFv4.3
+##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
+#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  M1  F1  C1
+chr1    100 snp03   A   G   .   .   .   GT  0/0 1/1 0/1
+chr1    200 snp05   C   T   .   .   .   GT  1/1 0/0 0/1
+chr1    300 snp09   G   A   .   .   .   GT  0/0 0/1 0/1
    +

    example.ped (six pedigree columns, of which only the +individual id in column 2 is used, then two alleles per marker)

    +
    FAM M1 0 0 1 -9 A A T T A A
+FAM F1 0 0 2 -9 G G C C A G
+FAM C1 M1 F1 0 -9 A G C T A G
    +

    example.map (one row per marker: chromosome, marker +name, genetic position, physical position)

    +
    1 snp03 0 100
+1 snp05 0 200
+1 snp09 0 300
    +

    Converting either file gives the same table:

    +
    geno_vcf = vcf_to_dosage("example.vcf", verbose = FALSE)   # ploidy read from the GT calls
+geno_vcf
    +
    + +
    +
    geno_ped = ped_to_dosage("example.ped", verbose = FALSE)   # example.map in the same folder gives the marker names
+geno_ped
    +
    + +
    +

    Both have the layout of genotypes and, here, the same +values:

    +
    genotypes[genotypes$id %in% c("M1", "F1", "C1"), c("id", "snp03", "snp05", "snp09")]
    +
    + +
    +

    In a .ped file the counted allele is the alphabetically +last allele at each marker (G, T and +G here). For breed composition, +format = "matrix" returns the same values with the ids as +row names instead of an id column:

    +
    geno_mat = vcf_to_dosage("example.vcf", format = "matrix", verbose = FALSE)
+geno_mat
    +
    ##    snp03 snp05 snp09
+## M1     0     2     0
+## F1     2     0     1
+## C1     1     1     1
+## attr(,"ploidy")
+## [1] 2
    +

    See ?vcf_to_dosage and ?ped_to_dosage for +the other options, such as reading a dosage field (DS), +call-rate filters and multiallelic markers.

    +
    +
    +

    3.2 Missing data

    +

    There are two kinds of missing information, coded differently:

    +
      +
    • Unknown parent (pedigree): 0. In a +fill_pedigree progeny table NA or an empty +cell also work.
      • +
    • Missing genotype call: NA in a dosage +table. The converters turn a VCF call containing . (for +example ./.) and a .ped allele coded +0 into NA.
      • +
    +

    Missing calls never count as errors. Each test only uses the markers +where every individual it compares has a call (the tests are explained +in How the Mendelian +error is calculated):

    +
      +
    • Trio test (validate_pedigree(), and +find_parentage() with best_pair or +fill_pedigree): markers where the offspring and both +parents have a call, the testable markers. A trio with fewer +testable markers than min_markers (default 10) is +low_markers; one with none is no_data.
      • +
    • Single-parent test (recommended corrections, and +the single-parent methods of find_parentage()): markers +where the parent and the offspring are both homozygous. If there are +none, the parent counts as failing (see C_low in Failed trios and +recommended correction).
      • +
    • Marker summary: mismatch_pct is +NA for a marker not tested in any trio.
      • +
    • Breed composition: a sample’s missing calls are +skipped marker by marker. A marker missing in every reference individual +of a line has no frequency (NaN in the output of +allele_freq_poly()), which makes +solve_composition_poly() stop with an error. Find these +markers with rownames(freq)[rowSums(is.na(freq)) > 0] +and remove them before rerunning.
      • +
    @@ -3761,102 +3849,295 @@

    4 Pedigree cleaning

    An individual is its own ancestor, or is recorded as its own parent Flagged only; must be fixed manually, since there is no safe -automatic choice +automatic choice. Each cycle is listed once, as a path in which every id +is a recorded parent of the one before it +(A -> B -> A) + + + +

    Unknown parents must be coded 0. With +verbose = TRUE the function prints a report of every +check:

    +
    clean_ped_results = check_ped(pedigree, verbose = TRUE)
    +
    ## 
+## === Pedigree Quality Check Report ===
+## 
+## Exact duplicate trios detected (removed in corrected pedigree):
+##    row_number        id male_parent female_parent
+## 10         10 C_missing          M1            F1
+## 11         11 C_missing          M1            F1
+## 
+## No conflicting trios found.
+## 
+## No missing parents found.
+## 
+## No IDs found as both male_parent and female_parent.
+## 
+## No dependencies detected.
+## 
+## The corrected pedigree is included in the returned list as corrected_pedigree.
    +

    The result is a list with one table per issue plus the corrected +pedigree. An empty table means the issue was not found:

    +
    names(clean_ped_results)
    +
    ## [1] "exact_duplicates"       "conflicting_trios"      "inconsistent_sex_roles"
+## [4] "missing_parents"        "dependencies"           "corrected_pedigree"
    +
    sapply(clean_ped_results, nrow)   # number of rows in each table
    +
    ##       exact_duplicates      conflicting_trios inconsistent_sex_roles 
+##                      2                      0                      0 
+##        missing_parents           dependencies     corrected_pedigree 
+##                      0                      0                     10
    + ++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
    ElementIn the toy data
    exact_duplicatesBoth copies of the repeated C_missing row; +row_number is their position in the input
    conflicting_triosEmpty: no id has two different sets of parents
    inconsistent_sex_rolesEmpty: no individual is used as both male and female parent
    missing_parentsEmpty: every parent has its own row
    dependenciesEmpty: no cycles
    corrected_pedigreeThe pedigree without the duplicate, ready for +validate_pedigree()
    -

    Unknown parents must be coded 0.

    -
    clean_ped_results = check_ped(pedigree, verbose = FALSE)
-
-# the corrected, analysis-ready pedigree
-clean_ped = clean_ped_results$corrected_pedigree
-clean_ped
    +
    clean_ped_results$exact_duplicates
    +
    + +
    +
    # the corrected, analysis-ready pedigree
+clean_ped = clean_ped_results$corrected_pedigree
+clean_ped
    -

    The result also holds one table per issue, for example -clean_ped_results$exact_duplicates (here the repeated -C_missing row) and -clean_ped_results$missing_parents. An empty table means the -issue was not found. Set verbose = TRUE to print a full -report to the console.

    5 Pedigree validation

    -

    validate_pedigree() compares every trio (offspring, male -parent, female parent) with the genotypes and computes a Mendelian error -rate. A marker is testable when the offspring and both parents -have a call.

    -
    ped_validate_results = validate_pedigree(clean_ped, genotypes,
-                                         verbose = FALSE, plot_results = TRUE)
    +

    validate_pedigree() tests every trio (offspring, male +parent, female parent) of the pedigree against the genotypes. It returns +the status of each trio, a recommended correction when a trio fails, a +corrected pedigree and a per-marker summary. Here the founders are +passed with founders_file, and +marker_trio_table = TRUE adds the marker-by-trio table used +below. Click the label under the code to see the console report.

    +
    ped_validate_results = validate_pedigree(clean_ped, genotypes,
+                                         founders_file     = founders,
+                                         marker_trio_table = TRUE,
+                                         verbose = TRUE, plot_results = TRUE)
    +
    + +Show the report printed with verbose = TRUE + +
    ## Found 1 trios with missing genotype data; flagged as no_genotype_data.
+## 
+## --- Trio Validation Summary ---
+## Total trios in pedigree: 10
+## fail                        : 1 (10.0%)
+## founders                    : 4 (40.0%)
+## low_markers                 : 1 (10.0%)
+## no_genotype_data            : 1 (10.0%)
+## pass                        : 3 (30.0%)
+## Error threshold: 5%
+## Homozygous threshold: 2%
+## Minimum markers required: 10
+## 
+## Correction summary:
+##   low_markers_remove_female_parent: 1
+##   remove_male_parent: 1
+## 
+##            id orig_male_parent orig_female_parent trio_mendelian_error_pct
+##        <char>           <char>             <char>                    <num>
+##  1:        M1                0                  0                       NA
+##  2:        F1                0                  0                       NA
+##  3:        M2                0                  0                       NA
+##  4:        F2                0                  0                       NA
+##  5:        C1               M1                 F1                     0.00
+##  6:        C2               M2                 F2                     0.00
+##  7:        C3               M1                 F2                     0.00
+##  8:     C_bad               M2                 F1                    41.67
+##  9:     C_low               M1                 F1                     0.00
+## 10: C_missing               M1                 F1                       NA
+##     trio_markers_tested           status           recommended_correction
+##                   <int>           <char>                           <char>
+##  1:                   0         founders                             none
+##  2:                   0         founders                             none
+##  3:                   0         founders                             none
+##  4:                   0         founders                             none
+##  5:                  12             pass                             none
+##  6:                  12             pass                             none
+##  7:                  12             pass                             none
+##  8:                  12             fail               remove_male_parent
+##  9:                   4      low_markers low_markers_remove_female_parent
+## 10:                   0 no_genotype_data                             none
+##     male_parent_hom_error_pct female_parent_hom_error_pct
+##                         <num>                       <num>
+##  1:                        NA                          NA
+##  2:                        NA                          NA
+##  3:                        NA                          NA
+##  4:                        NA                          NA
+##  5:                        NA                          NA
+##  6:                        NA                          NA
+##  7:                        NA                          NA
+##  8:                     83.33                           0
+##  9:                      0.00                          NA
+## 10:                        NA                          NA
    +

    -
    ped_validate_report = ped_validate_results$full_results
-ped_validate_report
    +

    The result is a list. full_results has one row per trio, +and pass, fail, low_markers, +no_genotype_data, founders and +missing_parents hold the same rows split by status. +corrected_pedigree, marker_summary and +marker_trio_table are described below, and +plot is the histogram above.

    +
    names(ped_validate_results)
    +
    ##  [1] "pass"               "fail"               "low_markers"       
+##  [4] "no_genotype_data"   "founders"           "missing_parents"   
+##  [7] "full_results"       "corrected_pedigree" "marker_summary"    
+## [10] "marker_trio_table"  "plot"
    +
    sapply(Filter(is.data.frame, ped_validate_results), nrow)   # rows in each table
    +
    ##               pass               fail        low_markers   no_genotype_data 
+##                  3                  1                  1                  1 
+##           founders    missing_parents       full_results corrected_pedigree 
+##                  4                  0                 10                 10 
+##     marker_summary  marker_trio_table 
+##                 12                 12
    +
    ped_validate_report = ped_validate_results$full_results
+ped_validate_report
    -
    -

    5.1 Parameters

    +

    In full_results, trio_mendelian_error_pct +is the trio error over trio_markers_tested testable +markers, and male_parent_hom_error_pct and +female_parent_hom_error_pct are the single-parent errors, +computed only for failed and low-marker trios. The thresholds used here +are the defaults; they and the other settings are arguments of the +function (see ?validate_pedigree).

    +
    +

    5.1 How the Mendelian +error is calculated

    +

    The trio error uses every marker where the offspring and both +parents have a call, heterozygous or homozygous, not only the +homozygous ones. At each testable marker the question is whether the two +parents can produce the offspring. Each parent passes on a gamete with +half of its alleles (ploidy / 2), so a parent with dosage g +passes between max(0, g − ploidy / 2) and min(g, +ploidy / 2) copies of allele B. The offspring dosage must lie between +the sum of the two parents’ minimums and the sum of their maximums; +otherwise the marker is a Mendelian error. For a diploid this gives the +familiar rules:

    ----- - - - + + + - - - + + + - - - + + + - - - + + + - - - + + + - - - + + + - - - - - - - - + + +
    ArgumentDefaultEffectMale parentFemale parentPossible offspring
    ploidy2Ploidy of the dosages000
    trio_error_threshold5Maximum % of testable markers with a Mendelian error for a trio to -pass010 or 1
    min_markers10Minimum testable markers needed to judge a trio; fewer gives -low_markers021
    single_parent_error_threshold2Maximum homozygous mismatch % for a single parent to pass the -single-parent check110, 1 or 2
    founders_fileNULLFile path, vector of ids or data.frame; ids listed here with both -parents 0 are reported as founders121 or 2
    marker_trio_tableFALSEAlso return the marker x trio table (see below)
    plot_resultsTRUEDraw a histogram of trio error percentages222
    +

    For a tetraploid, a parent with dosage 3 passes 1 or 2 copies, so a 3 +x 0 cross can only give 1 or 2. The trio error is then

    +
    +

    trio error (%) = 100 × markers with an error / testable markers

    +
    +

    and the trio passes when it is at or below +trio_error_threshold (default 5%). For odd ploidy +(e.g. triploid) gametes are unbalanced, so the trio test falls back to +the opposite-homozygote check of the single-parent test (see Failed trios and +recommended correction): a marker is an error when the offspring is +homozygous and a homozygous parent carries the other allele.

    +

    Here is the trio of C_bad (recorded as M2 x +F1) marker by marker, with the verdict from +marker_trio_table:

    +
    c_bad_trio = data.frame(
+  marker  = ped_validate_results$marker_trio_table$marker,
+  M2      = unlist(genotypes[genotypes$id == "M2",    -1]),
+  F1      = unlist(genotypes[genotypes$id == "F1",    -1]),
+  C_bad   = unlist(genotypes[genotypes$id == "C_bad", -1]),
+  verdict = ped_validate_results$marker_trio_table$C_bad_M2_F1,
+  row.names = NULL
+)
+c_bad_trio
    +
    + +
    +

    5 of the 12 testable markers are errors, a trio error of 41.67%, so +the trio fails. At snp01, for example, M2 is 2 +and F1 is 0, so every child must be 1, but +C_bad is 0. Heterozygous markers count too: at +snp10 the cross 2 x 1 can only give 1 or 2, and +C_bad is 0. All errors are labelled +mismatch_male, meaning M2 cannot be the parent +whatever F1 contributes; the single-parent test below +confirms this.

    5.2 Status of each @@ -3875,11 +4156,11 @@

    5.2 Status of each pass -Mendelian error at or below trio_error_threshold +Trio error at or below trio_error_threshold fail -Error above the threshold +Trio error above the threshold low_markers @@ -3888,37 +4169,57 @@

    5.2 Status of each no_genotype_data At least one of the three individuals is not in the genotype -file +table founders -Both parents 0 and the id is in +Both parents 0 and the id is listed in founders_file missing_both_parents, missing_male_parent, missing_female_parent -One or both parents are 0 +One or both parents are 0, and the individual is not a +listed founder no_data -No testable marker at all +No testable marker at all (only in full_results)

    In the toy data C1, C2 and C3 pass, C_bad fails, C_low is low_markers and C_missing has -no_genotype_data.

    +no_genotype_data. M1, F1, +M2 and F2 are founders because +they are listed in founders; without +founders_file they would be +missing_both_parents.

    5.4 Marker summary

    @@ -3985,7 +4296,7 @@

    5.4 Marker summary

    marker? A few markers with a mismatch % far above the rest (paralogous loci, poor clustering, null alleles) inflate the error of many trios. Remove them from the genotype file and rerun.

    -
    head(ped_validate_results$marker_summary)
    +
    head(ped_validate_results$marker_summary)
    -
    -
    -

    7.1 Parameters

    +

    When parentage is unknown, find_parentage() tests every +candidate in parents against each individual in +progeny and reports the best one. How candidates are scored +depends on the method; the score is reported as +mendelian_error_pct, and lower is better:

    @@ -4102,149 +4377,180 @@

    7.1 Parameters

    - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ArgumentDefaultEffect
    method"best_pair"The search (see below)
    ploidy2Ploidy of the dosages
    error_threshold5Error % above which the best candidate gets -high_error
    min_markers10If the best candidate has fewer testable markers, the status is -low_markers
    top_n1Report the best top_n candidates per progeny
    show_tiesTRUEReport tied_candidates
    allow_parent_selfingFALSEAllow the same individual as male and female parent
    exclude_self_matchTRUEDo not let a progeny match itself as parent
    founders_fileNULLFile path, vector of ids or data.frame; founders are skipped by -fill_pedigree
    - ---- - - - + + - - + + + - + + - + + - + +
    MethodSearchesCandidates testedScore
    best_pairMale x female pair with the lowest Mendelian error (needs both sexes -in the parents table)best_pair (default)Every male x female pairTrio error, as in validate_pedigree(): all markers +where the three individuals have a call
    best_male_parent, best_female_parentThe single candidate of that sex with the lowest homozygous mismatch -(compares only markers homozygous in both candidate and offspring); use -when the other parent is known or unavailableEach candidate of that sex aloneSingle-parent error: only markers homozygous in both candidate and +progeny
    best_matchSame homozygous mismatch test over candidates of either sexEach candidate of either sex aloneSingle-parent error, as above
    fill_pedigreeMissing parents of a 3-column pedigree (see below)The missing parent(s) of a pedigree; a known parent is kept +fixedTrio error (see Filling in missing +parents of a pedigree)
    -

    The best candidate is always reported. The status is -pass (error at or below error_threshold), -high_error (above it, so treat the assignment with -caution), low_markers, or -known_parent_no_genotype. In the toy data -C_bad, which failed validation, is recovered to its true -pair (M1/F1) at 0% error, while -C_low remains low_markers.

    +

    Only the single-parent methods rely on homozygous markers. Prefer +best_pair when both parents are among the candidates: it +uses all markers and both parents together, so it separates candidates +much better. Use the single-parent methods when the other parent is not +genotyped.

    +

    The best candidate is always reported, with a status: +pass (score at or below error_threshold, +default 5%), high_error (above it: no candidate fits well, +so the true parent may be missing from parents), +low_markers (fewer than min_markers markers +tested) or, for fill_pedigree, +known_parent_no_genotype. See ?find_parentage +for all options.

    +
    find_parentage_results = find_parentage(genotypes, parents, progeny,
+                                        method = "best_pair",
+                                        verbose = TRUE, plot_results = TRUE)
    +
    + +Show the report printed with verbose = TRUE + +
    ## Parent selfing is disallowed. Pairs with identical parents are removed.
+## 
+## === Parentage Assignment Report ===
+## 
+## Total progeny evaluated: 5
+## Method: best_pair  |  Error threshold: 5%  |  Minimum markers: 10
+## 
+## 4 progeny passed (80.0%).
+## 
+## No progeny flagged for high error.
+## 
+## 1 progeny flagged low_markers (20.0%):
+##        id male_parent female_parent mendelian_error_pct markers_tested
+##    <char>      <char>        <char>              <char>          <int>
+## 1:  C_low          M1            F1                0.00              4
+##         status tied_candidates
+##         <char>           <int>
+## 1: low_markers               1
+## 
+## Full results are included in the returned list as $full_results.
    +
    +

    +

    The result is a list like the one from +validate_pedigree(): full_results has one row +per progeny, pass, high_error and +low_markers hold the same rows split by status, and +plot is the histogram above.

    +
    names(find_parentage_results)
    +
    ## [1] "pass"         "high_error"   "low_markers"  "full_results" "plot"
    +
    sapply(Filter(is.data.frame, find_parentage_results), nrow)   # rows in each table
    +
    ##         pass   high_error  low_markers full_results 
+##            4            0            1            5
    +
    find_parentage_report = find_parentage_results$full_results
+find_parentage_report
    +
    +
    -
    -

    7.2 Top N and ties

    -

    tied_candidates is the number of candidates (pairs for -best_pair) that share the lowest error of a progeny; values -above 1 mean the best assignment is ambiguous. This is common with few -markers or with clones. With top_n > 1, rank 1 stays in -the base columns and ranks 2 to top_n are added as suffix -columns (male_parent_2, female_parent_2, -mendelian_error_pct_2).

    -
    find_parentage(genotypes, parents, progeny, method = "best_pair", top_n = 3,
-               verbose = FALSE, plot_results = FALSE)$full_results
    +

    Each row gives the best pair (male_parent, +female_parent), its error +(mendelian_error_pct), the number of markers tested, the +status and tied_candidates, the number of pairs that share +the best error. C1, C2 and C3 get +their recorded parents at 0% error, and C_bad, which failed +validation, is recovered to its true pair (M1 x +F1). C_low also gets M1 x +F1, but on only 4 markers, so it stays +low_markers. C_missing is not in the genotype +table, so it is dropped with a warning.

    +
    +

    6.1 Top N and ties

    +

    tied_candidates above 1 means several candidates share +the best error, so the assignment is ambiguous; this is common with few +markers or with clones. To see how close the alternatives are, set +top_n:

    +
    top3_results = find_parentage(genotypes, parents, progeny, method = "best_pair",
+                              top_n = 3, verbose = FALSE, plot_results = FALSE)
+names(top3_results)
    +
    ## [1] "pass"         "high_error"   "low_markers"  "full_results" "plot"
    +
    top3 = top3_results$full_results
+top3
    +

    The list has the same elements; only full_results gains +columns. Rank 1 stays in the base columns (male_parent, +female_parent, mendelian_error_pct, +markers_tested), and ranks 2 and 3 are added with the +suffixes _2 and _3 +(male_parent_2, female_parent_2, +mendelian_error_pct_2, markers_tested_2, and +so on). status and tied_candidates refer to +rank 1. A large gap between ranks means a clear assignment: for +C_bad the best pair has 0.00% error and the second best (M1 +x F2) 25%. With the default top_n = 1, tied best pairs are +added in the same suffix columns (show_ties = TRUE).

    -
    -

    7.3 Filling in missing +
    +

    6.2 Filling in missing parents of a pedigree

    With method = "fill_pedigree", progeny is a -pedigree with id, male_parent and -female_parent (for example -validate_pedigree()$corrected_pedigree). Parents coded -0, NA or empty are unknown. A known male -parent searches only females, a known female parent only males, and rows -with both unknown get the full pair search; rows with both parents known -are skipped. The parents table (with sex M, F -or A) is still required. The output adds -search_mode and known_parent_error_pct.

    -
    ped_to_fill = data.frame(id = c("C1", "C2"),
-                         male_parent = c("M1", "0"),
-                         female_parent = c("0", "0"))
-
-find_parentage(genotypes, parents, ped_to_fill, method = "fill_pedigree",
-               verbose = FALSE, plot_results = FALSE)$full_results
    +pedigree (id, male_parent, +female_parent), typically the +corrected_pedigree from validate_pedigree(). +Parents coded 0, NA or empty are unknown. Each +row is handled on its own: a known male parent searches only females, a +known female parent only males, and a row with both parents unknown gets +the full pair search, unless it is a listed founder. Rows with both +parents known are skipped. The parents table (with sex) is +still required.

    +
    filled_results = find_parentage(genotypes, parents,
+                                ped_validate_results$corrected_pedigree,
+                                method        = "fill_pedigree",
+                                founders_file = founders,
+                                verbose = FALSE, plot_results = FALSE)
+names(filled_results)
    +
    ## [1] "pass"         "high_error"   "low_markers"  "full_results" "plot"
    +
    filled_results$full_results
    +

    The list has the same elements as before. Only the rows with a parent +removed by validation are searched: the founders are skipped because +they are in founders, and the complete trios because both +parents are known. full_results gains two columns: +search_mode (known_male, +known_female, or pair when both parents are +unknown) and known_parent_error_pct, the single-parent +error of the parent that was kept (a high value suggests that it is +wrong too). C_bad kept F1, so only males were +tested, and M1 is found at 0%. C_low kept +M1 and gets F1 back, but on too few markers, +so it stays low_markers.

    -
    -

    8 Breed and line +
    +

    7 Breed and line composition

    BIGpopA estimates genome-wide breed or line composition for diploid and polyploid populations with a breedTools-style quadratic programming @@ -4253,26 +4559,26 @@

    8 Breed and line sample, and assign each sample to a line. This example mirrors an autotetraploid analysis (ploidy = 4) that assigns samples to two cultivar reference panels, Jewel and Draper.

    -
    # reference genotypes: individuals x SNPs, tetraploid dosage of allele B (0..4)
-reference = data.frame(
-  S1 = c(0, 1, 0, 4, 3, 4),
-  S2 = c(1, 0, 0, 3, 4, 4),
-  S3 = c(0, 1, 1, 4, 3, 3),
-  S4 = c(1, 0, 1, 3, 4, 3),
-  S5 = c(0, 1, 0, 4, 3, 4),
-  S6 = c(1, 0, 1, 3, 4, 3),
-  row.names = c("Jewel1","Jewel2","Jewel3","Draper1","Draper2","Draper3")
-)
-
-# reference panel membership (which individuals define each cultivar)
-ref_ids = list(
-  Jewel  = c("Jewel1","Jewel2","Jewel3"),
-  Draper = c("Draper1","Draper2","Draper3")
-)
-
-# allele frequencies of the reference panels (SNPs x panels)
-freq = allele_freq_poly(reference, ref_ids, ploidy = 4)
-freq
    +
    # reference genotypes: individuals x SNPs, tetraploid dosage of allele B (0..4)
+reference = data.frame(
+  S1 = c(0, 1, 0, 4, 3, 4),
+  S2 = c(1, 0, 0, 3, 4, 4),
+  S3 = c(0, 1, 1, 4, 3, 3),
+  S4 = c(1, 0, 1, 3, 4, 3),
+  S5 = c(0, 1, 0, 4, 3, 4),
+  S6 = c(1, 0, 1, 3, 4, 3),
+  row.names = c("Jewel1","Jewel2","Jewel3","Draper1","Draper2","Draper3")
+)
+
+# reference panel membership (which individuals define each cultivar)
+ref_ids = list(
+  Jewel  = c("Jewel1","Jewel2","Jewel3"),
+  Draper = c("Draper1","Draper2","Draper3")
+)
+
+# allele frequencies of the reference panels (SNPs x panels)
+freq = allele_freq_poly(reference, ref_ids, ploidy = 4)
+freq
    ##         Jewel    Draper
 ## S1 0.08333333 0.9166667
 ## S2 0.08333333 0.9166667
@@ -4283,28 +4589,39 @@ 
    8 Breed and line
 
    The references should be representative, unadmixed individuals; the
 more distinct the references, the better the estimate. Marker names must
 match between the reference and the samples.
    
-
    # samples to assign: samples x SNPs (same SNPs as the reference)
-validation = data.frame(
-  S1 = c(0, 4, 2),
-  S2 = c(1, 3, 2),
-  S3 = c(0, 4, 2),
-  S4 = c(1, 3, 2),
-  S5 = c(0, 4, 2),
-  S6 = c(1, 3, 2),
-  row.names = c("Sample1","Sample2","Sample3")
-)
-
-# composition (one column per panel plus an R2 model-fit column)
-prediction = as.data.frame(solve_composition_poly(validation, freq, ploidy = 4))
-prediction
    
+
    # samples to assign: samples x SNPs (same SNPs as the reference)
+validation = data.frame(
+  S1 = c(0, 4, 2),
+  S2 = c(1, 3, 2),
+  S3 = c(0, 4, 2),
+  S4 = c(1, 3, 2),
+  S5 = c(0, 4, 2),
+  S6 = c(1, 3, 2),
+  row.names = c("Sample1","Sample2","Sample3")
+)
+
+# composition (one column per panel plus an R2 model-fit column)
+prediction = as.data.frame(solve_composition_poly(validation, freq, ploidy = 4))
+prediction
    
 
    
 
 
    
-
    The model returns an R2 value, but it is expected to be low because a
-discrete allele dosage is regressed on a continuous allele
-frequency.
    
+
    R2 measures how well the estimated mixture of reference
+lines explains a sample’s genotypes. The sample’s dosages are divided by
+the ploidy (0, 0.25, …, 1 for a tetraploid); the predicted value at each
+marker is the allele frequency of the reference lines weighted by the
+estimated proportions; and R2 is the squared correlation
+between the two across the sample’s non-missing markers. It describes
+the fit, not the accuracy of the proportions. It is usually well below 1
+even for a pure individual, because each genotype is one of a few
+discrete values while the prediction is a frequency, so every marker at
+which the line is polymorphic adds unexplained variation. Compare
+R2 between samples analysed with the same reference panel
+and markers: a sample with a much lower R2 than the rest
+may come from a population that is not in the reference panel, or have
+poor-quality genotypes.
    
 
    assign_composition_line() turns the result into the
 table used by Familia: the estimated percentage of each panel and a
 Predicted line. A sample is assigned to its
@@ -4312,8 +4629,8 @@ 
    8 Breed and line
 assign_threshold (default 50); otherwise it is labelled
 Undetermined. With assign_threshold = 0 every
 sample is assigned.
    
-
    pred_results = assign_composition_line(solve_composition_poly(validation, freq, ploidy = 4))
-pred_results
    
+
    pred_results = assign_composition_line(solve_composition_poly(validation, freq, ploidy = 4))
+pred_results