From da4162e8177e8f2de4e1583a2da1e603cd12f75f Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Mon, 5 Oct 2026 15:32:09 -0400 Subject: [PATCH 1/8] fix for missrecorded selfing --- DESCRIPTION | 2 +- NEWS.md | 4 +++ R/validate_pedigree.R | 35 ++++++++++++++++++++++--- tests/testthat/test-validate_pedigree.R | 22 +++++++++++++++- 4 files changed, 58 insertions(+), 5 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 5e9c25b..2c9f235 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BIGpopA Title: Pedigree Validation Genetic Composition of Diploids & Polyploids -Version: 2.1.0 +Version: 2.1.1 Authors@R: c(person(given='Josue', family='Chinchilla-Vargas', email='josue.chinchilla@ufl.edu', diff --git a/NEWS.md b/NEWS.md index f8366fc..3bbbe8f 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,7 @@ +# BIGpopA 2.1.1 +- `validate_pedigree()`: a trio that fails the Mendelian error threshold while neither parent fails the single-parent homozygous check now gets `recommended_correction = "unresolved"` instead of `"keep_both"`. This happens, for example, with a recorded self whose offspring carries alleles the parent lacks: the pedigree is wrong, but the faulty parent cannot be identified. These trios are left unchanged in `corrected_pedigree` for manual review. +- `validate_pedigree()` plot: bars are coloured by trio error first, so no trio above the threshold is shown as Pass. Unresolved trios are shown in grey ("Fail - Parents Unresolved"), and the legend wraps onto two rows so all labels are visible. + # BIGpopA 2.1.0 - New `vcf_to_dosage()` converts a VCF (`.vcf` / `.vcf.gz`) or `vcfR` object into allele-B dosages for any ploidy, as an `id` + markers table or an individuals x markers matrix. - New `ped_to_dosage()` converts a PLINK `.ped` file (diploid) into allele dosages. A `.map` file is optional and only supplies marker names. The counted allele is the alphabetically/numerically last allele at each marker, or can be supplied to code a second file like a first one. diff --git a/R/validate_pedigree.R b/R/validate_pedigree.R index d0e7046..38a7c12 100644 --- a/R/validate_pedigree.R +++ b/R/validate_pedigree.R @@ -50,6 +50,26 @@ #' \item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} #' } #' +#' @details The `recommended_correction` column of `full_results` takes these +#' values: +#' \describe{ +#' \item{none}{Trio passes, or a parent is missing / no genotype data.} +#' \item{remove_male_parent, remove_female_parent}{The trio fails and only one +#' parent fails the single-parent homozygous check; that parent is set to 0 +#' in `corrected_pedigree`.} +#' \item{remove_both}{The trio fails and both parents fail the homozygous +#' check; both are set to 0.} +#' \item{unresolved}{The trio fails but neither parent fails the homozygous +#' check (for example a recorded self whose offspring carries alleles the +#' parent lacks). The pedigree is likely wrong but the faulty parent cannot +#' be identified, so `corrected_pedigree` is left unchanged for manual review.} +#' \item{keep_both}{Only with low markers: the trio is within the threshold +#' and both parents pass.} +#' } +#' Low-marker trios get the same values prefixed with `low_markers_`. In the +#' plot, every trio above the error threshold is coloured as a failure; +#' unresolved trios are shown in grey. +#' #' @examples #' \donttest{ #' geno_df <- data.frame( @@ -287,7 +307,11 @@ validate_pedigree <- function(pedigree_file, genotypes_file, female_acceptable <- !is.na(female_parent_error_pct) && female_parent_error_pct <= single_parent_error_threshold if (male_acceptable && female_acceptable) { - correction_decision <- "keep_both" + # Neither parent is excluded by the homozygous check. If the trio + # still fails (e.g. a recorded self whose offspring carries alleles + # absent from the parent), the pedigree is wrong but the faulty + # parent cannot be identified: flag it instead of keeping it. + correction_decision <- if (error_pct > trio_error_threshold) "unresolved" else "keep_both" } else if (male_acceptable && !female_acceptable) { correction_decision <- "remove_female_parent" best_f <- find_best_parent(prog_id, exclude_ids = c(male_parent_id)) @@ -407,15 +431,17 @@ validate_pedigree <- function(pedigree_file, genotypes_file, warning("ggplot2 is required for plot_results = TRUE. Please install it.", call. = FALSE) } else { plot_df <- final_df[!is.na(final_df$trio_mendelian_error_pct)] + # Colour by the trio error first, so no trio above the threshold is shown as pass plot_df$plot_status <- dplyr::case_when( - plot_df$recommended_correction %in% c("none", "keep_both", - "low_markers_keep_both") ~ "pass", + plot_df$trio_mendelian_error_pct <= trio_error_threshold ~ "pass", plot_df$recommended_correction %in% c("remove_male_parent", "remove_female_parent", "low_markers_remove_male_parent", "low_markers_remove_female_parent") ~ "fail_one_parent", plot_df$recommended_correction %in% c("remove_both", "low_markers_remove_both") ~ "fail_both_parents", + plot_df$recommended_correction %in% c("unresolved", + "low_markers_unresolved") ~ "fail_unresolved", TRUE ~ "other" ) n_total <- nrow(plot_df) @@ -438,12 +464,15 @@ validate_pedigree <- function(pedigree_file, genotypes_file, values = c("pass" = "#339900", "fail_one_parent" = "#F1C40F", "fail_both_parents" = "#cc3333", + "fail_unresolved" = "#808080", "other" = "#BDC3C7"), labels = c("pass" = "Pass", "fail_one_parent" = "Fail - One Parent", "fail_both_parents" = "Fail - Both Parents", + "fail_unresolved" = "Fail - Parents Unresolved", "other" = "Other") ) + + ggplot2::guides(fill = ggplot2::guide_legend(nrow = 2)) + # keeps every label visible ggplot2::labs( title = "Trio Mendelian Error Distribution", subtitle = paste0("Trios with Genotype Data Tested: ", n_total, diff --git a/tests/testthat/test-validate_pedigree.R b/tests/testthat/test-validate_pedigree.R index 7930ae2..520bd62 100644 --- a/tests/testthat/test-validate_pedigree.R +++ b/tests/testthat/test-validate_pedigree.R @@ -262,7 +262,27 @@ test_that("fail trio with one acceptable parent gets remove_* correction", { r <- out$full_results[id == "IND_D"] expect_true(r$recommended_correction %in% c("remove_male_parent", "remove_female_parent", "remove_both", - "keep_both")) + "unresolved")) +}) + +test_that("a failing trio whose parents both pass the homozygous check is unresolved", { + # Recorded self: P1 x P1, but the offspring is heterozygous where P1 is + # homozygous at 10 of 50 markers, so it cannot be a self of P1. + p1 <- base::rep(c(0L, 2L), 25) + off <- p1 + off[1:10] <- 1L + geno <- base::data.frame(id = c("P1", "Off"), base::rbind(p1, off), row.names = NULL) + base::names(geno)[-1] <- base::paste0("SNP", 1:50) + ped <- base::data.frame(id = "Off", male_parent = "P1", female_parent = "P1") + + out <- validate_pedigree(ped, geno, verbose = FALSE, plot_results = TRUE) + r <- out$full_results[id == "Off"] + expect_equal(r$status, "fail") + expect_equal(r$recommended_correction, "unresolved") + # Pedigree left unchanged for manual review + expect_equal(out$corrected_pedigree[id == "Off"]$male_parent, "P1") + # Plotted as a failure (grey), not as pass + expect_equal(out$plot$data$plot_status[out$plot$data$id == "Off"], "fail_unresolved") }) test_that("trio_mendelian_error_pct is 0 for a perfect Mendelian trio", { From 59fca99488bbbad26393420894756c5ff7a05e71 Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Mon, 5 Oct 2026 15:32:41 -0400 Subject: [PATCH 2/8] updated manual --- man/validate_pedigree.Rd | 21 +++++++++++++++++++++ 1 file changed, 21 insertions(+) diff --git a/man/validate_pedigree.Rd b/man/validate_pedigree.Rd index 0954eae..85655f1 100644 --- a/man/validate_pedigree.Rd +++ b/man/validate_pedigree.Rd @@ -77,6 +77,27 @@ best-matching replacements, and outputs a corrected pedigree. Founder trios supplied. Trios absent from the genotype file are retained as no_genotype_data. } +\details{ +The \code{recommended_correction} column of \code{full_results} takes these +values: +\describe{ +\item{none}{Trio passes, or a parent is missing / no genotype data.} +\item{remove_male_parent, remove_female_parent}{The trio fails and only one +parent fails the single-parent homozygous check; that parent is set to 0 +in \code{corrected_pedigree}.} +\item{remove_both}{The trio fails and both parents fail the homozygous +check; both are set to 0.} +\item{unresolved}{The trio fails but neither parent fails the homozygous +check (for example a recorded self whose offspring carries alleles the +parent lacks). The pedigree is likely wrong but the faulty parent cannot +be identified, so \code{corrected_pedigree} is left unchanged for manual review.} +\item{keep_both}{Only with low markers: the trio is within the threshold +and both parents pass.} +} +Low-marker trios get the same values prefixed with \code{low_markers_}. In the +plot, every trio above the error threshold is coloured as a failure; +unresolved trios are shown in grey. +} \examples{ \donttest{ geno_df <- data.frame( From b7d2ee3c9ca77d462de5357cbe6451865473b23c Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Wed, 7 Oct 2026 11:44:26 -0400 Subject: [PATCH 3/8] find_parentage(): New top_n argument (best N candidates per progeny) and a tied_candidates column that flags ambiguous assignments. New method = "fill_pedigree" fills missing parents of a 3-column pedigree, with a founders_file option. It errors if the input isn't such a pedigree. validate_pedigree(): New marker_summary with per-marker mismatch counts. Suggested replacement parents (best_*_candidate) removed. recommended_correction and corrected_pedigree unchanged. Tests: oracle-based tests for find_parentage() passed (957). The marker_summary tests passed (25). Not re-run since the suggested-parent removal. Docs: .Rd, vignette, README, NEWS.md, cran-comments.md and DESCRIPTION (2.2.0) updated. --- DESCRIPTION | 2 +- NEWS.md | 7 + R/find_parentage.R | 253 ++++++++++++++++-- R/validate_pedigree.R | 102 +++---- README.md | 4 +- cran-comments.md | 31 +-- man/find_parentage.Rd | 49 +++- man/validate_pedigree.Rd | 9 +- tests/testthat/test-find_parentage.R | 148 ++++++++++ tests/testthat/test-validate_pedigree.R | 34 +-- .../test-validate_pedigree_marker_summary.R | 113 ++++++++ vignettes/BIGpopA.Rmd | 50 +++- 12 files changed, 643 insertions(+), 159 deletions(-) create mode 100644 tests/testthat/test-validate_pedigree_marker_summary.R diff --git a/DESCRIPTION b/DESCRIPTION index 2c9f235..c55da55 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: BIGpopA Title: Pedigree Validation Genetic Composition of Diploids & Polyploids -Version: 2.1.1 +Version: 2.2.0 Authors@R: c(person(given='Josue', family='Chinchilla-Vargas', email='josue.chinchilla@ufl.edu', diff --git a/NEWS.md b/NEWS.md index 3bbbe8f..86483ce 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,10 @@ +# BIGpopA 2.2.0 +- `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. +- `find_parentage()`: new `method = "fill_pedigree"` fills in missing parents of a pedigree. `progeny_file` must have `id`, `male_parent` and `female_parent` columns (e.g. `validate_pedigree()$corrected_pedigree`); parents coded `0`/`NA`/empty are unknown. A known male searches only females, a known female only males, and rows with both unknown get the full pair search (skip founders with the new `founders_file`). Rows with both parents known are skipped. Results gain `search_mode` and `known_parent_error_pct`, and a known parent without genotypes gets `status = "known_parent_no_genotype"`. +- `find_parentage()`: new `top_n` argument reports the best `top_n` candidates per progeny, ranked by Mendelian error (ties broken by markers tested). 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`, `markers_tested_2` for `best_pair`; `best_match_2` and the same error/marker columns for the other methods). The default `top_n = 1` keeps the previous output, including `show_ties`. +- `find_parentage()`: new `tied_candidates` column gives the number of candidates (pairs for `best_pair`) that share each progeny's lowest error. Values above 1 mean the best assignment is ambiguous, including for the single-parent methods, which previously did not report ties. + # BIGpopA 2.1.1 - `validate_pedigree()`: a trio that fails the Mendelian error threshold while neither parent fails the single-parent homozygous check now gets `recommended_correction = "unresolved"` instead of `"keep_both"`. This happens, for example, with a recorded self whose offspring carries alleles the parent lacks: the pedigree is wrong, but the faulty parent cannot be identified. These trios are left unchanged in `corrected_pedigree` for manual review. - `validate_pedigree()` plot: bars are coloured by trio error first, so no trio above the threshold is shown as Pass. Unresolved trios are shown in grey ("Fail - Parents Unresolved"), and the legend wraps onto two rows so all labels are visible. diff --git a/R/find_parentage.R b/R/find_parentage.R index 4d1ab0f..e679ff9 100644 --- a/R/find_parentage.R +++ b/R/find_parentage.R @@ -14,9 +14,12 @@ #' data.table with an 'id' column and an optional 'sex' column #' ('M', 'F', or 'A'). If absent, all parents are treated as ambiguous. #' @param progeny_file Path to a TSV/CSV/TXT file, OR a data.frame / -#' data.table with an 'id' column. +#' data.table with an 'id' column. For `method = "fill_pedigree"` it must +#' be a pedigree with `id`, `male_parent` and `female_parent` columns (e.g. +#' the `corrected_pedigree` from [validate_pedigree()]); see Details. #' @param method Character. One of "best_male_parent", "best_female_parent", -#' "best_match", or "best_pair" (default). +#' "best_match", "best_pair" (default), or "fill_pedigree" (requires a +#' pedigree as `progeny_file`; see Details). #' @param min_markers Integer. Minimum markers required; fewer flags #' low_markers (default: 10). #' @param error_threshold Numeric. Maximum mismatch percentage; exceeded values @@ -39,13 +42,46 @@ #' Odd ploidy (e.g. triploid), where balanced gametes are undefined, falls #' back to a model-free opposite-homozygote exclusion evaluated on #' homozygous-informative markers only (reduced power). Default is 2. +#' @param founders_file Character, optional. Path to a one-column file of +#' founder IDs. Only used when `progeny_file` is a pedigree: founders with +#' both parents unknown are skipped instead of searched. Default is NULL. +#' @param top_n Integer >= 1. Number of best candidates to report per progeny. +#' With the default of 1, only the best assignment is reported (plus tied +#' best pairs when `show_ties = TRUE` and `method = "best_pair"`). With +#' `top_n > 1`, candidates are ranked by error (ties broken by the number of +#' markers tested) and the best `top_n` are reported: rank 1 in the base +#' columns and ranks 2 to `top_n` in suffix columns (`male_parent_2`, +#' `female_parent_2`, `mendelian_error_pct_2`, `markers_tested_2` for +#' "best_pair"; `best_match_2`, `mendelian_error_pct_2`, `markers_tested_2` +#' for the other methods). Ties are then ranked like any other candidate, so +#' `show_ties` has no effect. Progeny with fewer than `top_n` candidates +#' have `NA` in the remaining ranks. Status refers to rank 1 only. +#' +#' @details **Method "fill_pedigree" ("fill in the blanks").** The +#' `progeny_file` must have `male_parent` and `female_parent` columns. A parent coded `0`, `NA` or +#' empty is treated as unknown and each row is handled on its own: +#' a known male searches only females, a known female searches only males, +#' and a row with both unknown gets the full pair search (unless listed in +#' `founders_file`). Rows with both parents known are skipped. The search is +#' the "best_pair" Mendelian trio test with the known parent held +#' fixed, and candidates come from `parents_file` (sex M/F/A as before); the +#' known parent only needs genotypes. Two columns are added to the results: +#' `search_mode` ("pair", "known_male" or "known_female") and +#' `known_parent_error_pct` (homozygous mismatch between the known parent and +#' the progeny; a high value suggests the known parent is wrong). A known +#' parent absent from the genotype file gives +#' `status = "known_parent_no_genotype"`. Without `method = "fill_pedigree"`, +#' any `male_parent` / `female_parent` columns in `progeny_file` are ignored. #' #' @return A named list (returned invisibly) with elements: #' \describe{ #' \item{pass}{Progeny with a confident parentage assignment.} #' \item{high_error}{Progeny whose best assignment exceeds the error threshold.} #' \item{low_markers}{Progeny with insufficient markers for a valid assignment.} -#' \item{full_results}{Complete data.table with all progeny and all output columns.} +#' \item{full_results}{Complete data.table with all progeny and all output columns. +#' `tied_candidates` gives the number of candidates (pairs for "best_pair") +#' that share the lowest error for each progeny; values above 1 mean the +#' best assignment is ambiguous.} #' \item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} #' } #' @@ -101,7 +137,9 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, exclude_self_match = TRUE, verbose = TRUE, plot_results = TRUE, - ploidy = 2) { + ploidy = 2, + top_n = 1, + founders_file = NULL) { ## silence R CMD check NOTEs id <- sex <- male_parent <- female_parent <- NULL @@ -109,7 +147,7 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, #### Input Validation #### allowed_methods <- c("best_male_parent", "best_female_parent", - "best_match", "best_pair") + "best_match", "best_pair", "fill_pedigree") if (!method %in% allowed_methods) stop("Method must be one of: ", paste(allowed_methods, collapse = ", ")) if (min_markers < 1) @@ -117,6 +155,10 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, if (error_threshold < 0 || error_threshold > 100) stop("error_threshold must be between 0 and 100.") .check_ploidy(ploidy) + if (!is.numeric(top_n) || length(top_n) != 1 || is.na(top_n) || + top_n < 1 || top_n != base::round(top_n)) + stop("top_n must be a single positive integer.") + top_n <- base::as.integer(top_n) # Accept file path OR in-memory data.frame / data.table read_flex <- function(x, label, ...) { @@ -143,7 +185,35 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, # Column names are matched ignoring case (e.g. ID, Sex) all_parents <- .standardize_names(all_parents, c("id", "sex")) - progeny_candidates <- .standardize_names(progeny_candidates, "id") + progeny_candidates <- .standardize_names(progeny_candidates, + c("id", "male_parent", "female_parent")) + + #### Pedigree ("fill in the blanks") mode #### + is_unknown <- function(x) base::is.na(x) | base::trimws(base::as.character(x)) %in% c("", "0", "NA") + method_reported <- method + pedigree_mode <- method == "fill_pedigree" + if (pedigree_mode) { + if (!base::all(c("male_parent", "female_parent") %in% base::names(progeny_candidates))) + stop("method = 'fill_pedigree' requires a progeny_file with id, male_parent ", + "and female_parent columns.") + method <- "best_pair" + progeny_candidates[, male_parent := base::as.character(male_parent)] + progeny_candidates[, female_parent := base::as.character(female_parent)] + founder_ids <- if (!is.null(founders_file)) + base::unique(data.table::fread(founders_file, header = FALSE, + colClasses = "character")[[1]]) else + base::character(0) + um <- is_unknown(progeny_candidates$male_parent) + uf <- is_unknown(progeny_candidates$female_parent) + is_founder <- um & uf & progeny_candidates$id %in% founder_ids + if (verbose) + base::cat("Pedigree mode:", base::sum(!um & !uf), "complete trios skipped,", + base::sum(is_founder), "founders skipped,", + base::sum((um | uf) & !is_founder), "rows with a missing parent to search.\n") + progeny_candidates <- progeny_candidates[(um | uf) & !is_founder] + if (base::nrow(progeny_candidates) == 0) + stop("No progeny with a missing parent to fill in.") + } valid_ids <- genos$id removed_parents <- base::setdiff(all_parents$id, valid_ids) @@ -211,8 +281,10 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, best_match = NA_character_, mendelian_error_pct = NA_real_, markers_tested = NA_integer_, - status = NA_character_ + status = NA_character_, + tied_candidates = NA_integer_ ) + extra_list <- base::vector("list", n_progeny) for (i in base::seq_len(n_progeny)) { prog_id <- base::rownames(progeny_genos)[i] @@ -231,7 +303,15 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, percent_mismatch[self_idx] <- NA_real_ } - best_idx <- base::which.min(percent_mismatch) + if (top_n > 1) { + valid_idx <- base::which(!base::is.na(percent_mismatch)) + ranked_idx <- valid_idx[base::order(percent_mismatch[valid_idx], + -comparisons[valid_idx])] + top_idx <- ranked_idx[base::seq_len(base::min(top_n, base::length(ranked_idx)))] + best_idx <- top_idx[base::seq_len(base::min(1L, base::length(top_idx)))] + } else { + best_idx <- base::which.min(percent_mismatch) + } if (base::length(best_idx) == 0) { data.table::set(results_dt, i, "markers_tested", 0L) data.table::set(results_dt, i, "status", "low_markers") @@ -244,6 +324,35 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, data.table::set(results_dt, i, "mendelian_error_pct", best_error) data.table::set(results_dt, i, "markers_tested", base::as.integer(best_markers)) data.table::set(results_dt, i, "status", assign_status(best_markers, best_error)) + data.table::set(results_dt, i, "tied_candidates", + base::as.integer(base::sum(percent_mismatch == percent_mismatch[best_idx], + na.rm = TRUE))) + + if (top_n > 1 && base::length(top_idx) > 1) { + rank_row <- base::list(id = prog_id) + for (k in base::seq(2, base::length(top_idx))) { + rank_row[[base::paste0("best_match_", k)]] <- base::rownames(parent_genos)[top_idx[k]] + rank_row[[base::paste0("mendelian_error_pct_", k)]] <- base::round(percent_mismatch[top_idx[k]], 2) + rank_row[[base::paste0("markers_tested_", k)]] <- base::as.integer(comparisons[top_idx[k]]) + } + extra_list[[i]] <- data.table::as.data.table(rank_row) + } + } + + # Add ranks 2..top_n as suffix columns, keeping the progeny order + extra_dt <- data.table::rbindlist(base::Filter(Negate(base::is.null), extra_list), + fill = TRUE, use.names = TRUE) + if (top_n > 1) { + for (k in base::seq(2, top_n)) { + cols <- base::paste0(c("best_match_", "mendelian_error_pct_", "markers_tested_"), k) + for (col in cols) { + vals <- if (base::nrow(extra_dt) > 0 && col %in% base::names(extra_dt)) + extra_dt[[col]][base::match(results_dt$id, extra_dt$id)] else + base::rep(if (grepl("^best_match", col)) NA_character_ else + if (grepl("^markers", col)) NA_integer_ else NA_real_, base::nrow(results_dt)) + data.table::set(results_dt, j = col, value = vals) + } + } } final_df <- results_dt } @@ -261,7 +370,20 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, male_parent_genos_mat <- genos_mat[parent_pairs$male_parent, , drop = FALSE] female_parent_genos_mat <- genos_mat[parent_pairs$female_parent, , drop = FALSE] progeny_ids <- progeny_candidates$id - progeny_mat <- genos_mat[progeny_ids, , drop = FALSE] + n_all <- base::length(progeny_ids) + if (pedigree_mode) { + kmale <- base::ifelse(is_unknown(progeny_candidates$male_parent), + NA_character_, progeny_candidates$male_parent) + kfem <- base::ifelse(is_unknown(progeny_candidates$female_parent), + NA_character_, progeny_candidates$female_parent) + } else { + kmale <- kfem <- base::rep(NA_character_, n_all) + } + search_mode <- base::ifelse(!base::is.na(kmale), "known_male", + base::ifelse(!base::is.na(kfem), "known_female", "pair")) + pair_idx <- base::which(search_mode == "pair") + # Full pair matrices are only built for progeny with no known parent + progeny_mat <- genos_mat[progeny_ids[pair_idx], , drop = FALSE] n_progeny <- base::nrow(progeny_mat) n_pairs <- base::nrow(parent_pairs) @@ -305,18 +427,80 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, female_parent = NA_character_, mendelian_error_pct = NA_character_, markers_tested = NA_integer_, - status = NA_character_ + status = NA_character_, + tied_candidates = NA_integer_ ) - results_list <- base::vector("list", n_progeny) - for (j in base::seq_len(n_progeny)) { - prog_id <- progeny_ids[j] - percent_mismatch <- pct_mismatch_mat[, j] - comparisons <- comparison_mat[, j] + # Known-parent bookkeeping (pedigree mode only) + known_id_vec <- base::ifelse(search_mode == "known_male", kmale, + base::ifelse(search_mode == "known_female", kfem, + NA_character_)) + if (pedigree_mode) { + known_err <- base::rep(NA_real_, n_all) + for (j in base::which(search_mode != "pair")) { + kid <- known_id_vec[j] + if (kid %in% base::rownames(genos_hom_mat)) { + kh <- genos_hom_mat[kid, ] + ph <- genos_hom_mat[progeny_ids[j], ] + n_cmp <- base::sum(!base::is.na(kh) & !base::is.na(ph)) + known_err[j] <- if (n_cmp == 0) NA_real_ else + base::round(base::sum(kh != ph, na.rm = TRUE) / n_cmp * 100, 2) + } + } + data.table::set(results_dt, j = "search_mode", value = search_mode) + data.table::set(results_dt, j = "known_parent_error_pct", value = known_err) + data.table::set(results_dt, which(search_mode == "known_male"), + "male_parent", known_id_vec[search_mode == "known_male"]) + data.table::set(results_dt, which(search_mode == "known_female"), + "female_parent", known_id_vec[search_mode == "known_female"]) + } + + results_list <- base::vector("list", n_all) + for (j in base::seq_len(n_all)) { + prog_id <- progeny_ids[j] + mode_j <- search_mode[j] + + if (mode_j == "pair") { + pairs_j <- parent_pairs + pj <- base::match(j, pair_idx) + percent_mismatch <- pct_mismatch_mat[, pj] + comparisons <- comparison_mat[, pj] + } else { + known_id <- known_id_vec[j] + if (!known_id %in% base::rownames(genos_mat)) { + warning("Known parent '", known_id, "' of progeny '", prog_id, + "' is not in the genotype file; this progeny cannot be tested.", + call. = FALSE) + data.table::set(results_dt, j, "markers_tested", 0L) + data.table::set(results_dt, j, "status", "known_parent_no_genotype") + next + } + cand <- if (mode_j == "known_male") female_parent_candidates$id else + male_parent_candidates$id + if (!allow_parent_selfing) cand <- cand[cand != known_id] + if (exclude_self_match) cand <- cand[cand != prog_id] + cand <- cand[cand %in% base::rownames(genos_mat)] + if (base::length(cand) == 0) { + data.table::set(results_dt, j, "markers_tested", 0L) + data.table::set(results_dt, j, "status", "low_markers") + next + } + pairs_j <- if (mode_j == "known_male") + data.table::data.table(male_parent = known_id, female_parent = cand) else + data.table::data.table(male_parent = cand, female_parent = known_id) + pm <- genos_mat[pairs_j$male_parent, , drop = FALSE] + pf <- genos_mat[pairs_j$female_parent, , drop = FALSE] + pg <- base::matrix(genos_mat[prog_id, ], nrow = base::nrow(pairs_j), + ncol = base::ncol(genos_mat), byrow = TRUE) + mm <- base::rowSums(.mend_mismatch(pm, pf, pg, ploidy), na.rm = TRUE) + comparisons <- base::rowSums(.mend_testable(pm, pf, pg, ploidy)) + percent_mismatch <- (mm / comparisons) * 100 + percent_mismatch[base::is.nan(percent_mismatch)] <- NA + } if (exclude_self_match) { - self_pair_idx <- parent_pairs$male_parent == prog_id | - parent_pairs$female_parent == prog_id + self_pair_idx <- pairs_j$male_parent == prog_id | + pairs_j$female_parent == prog_id percent_mismatch[self_pair_idx] <- NA_real_ } @@ -328,39 +512,52 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, next } - best_indices <- base::which(percent_mismatch == min_mismatch_val) - if (base::length(best_indices) > 1) { - best_markers_per_pair <- comparisons[best_indices] - max_markers <- base::max(best_markers_per_pair) - best_indices <- best_indices[best_markers_per_pair == max_markers] + if (top_n > 1) { + # Rank all candidate pairs by error, then by markers tested; keep top_n + valid_idx <- base::which(!base::is.na(percent_mismatch)) + ranked_idx <- valid_idx[base::order(percent_mismatch[valid_idx], + -comparisons[valid_idx])] + best_indices <- ranked_idx[base::seq_len(base::min(top_n, base::length(ranked_idx)))] + } else { + best_indices <- base::which(percent_mismatch == min_mismatch_val) + if (base::length(best_indices) > 1) { + best_markers_per_pair <- comparisons[best_indices] + max_markers <- base::max(best_markers_per_pair) + best_indices <- best_indices[best_markers_per_pair == max_markers] + } } - best_pairs <- parent_pairs[best_indices] + best_pairs <- pairs_j[best_indices] best_markers <- comparisons[best_indices[1]] best_error <- base::round(min_mismatch_val, 2) a_status <- assign_status(best_markers, best_error) - if (!show_ties && base::nrow(best_pairs) > 1) { + if (top_n == 1 && !show_ties && base::nrow(best_pairs) > 1) { warning("Progeny '", prog_id, "' has ", base::nrow(best_pairs), " tied best pairs. Only one is reported as show_ties=FALSE.", call. = FALSE) } + report_all <- show_ties || top_n > 1 num_to_report <- base::min(base::nrow(best_pairs), - if (show_ties) base::nrow(best_pairs) else 1) + if (report_all) base::nrow(best_pairs) else 1) data.table::set(results_dt, j, "male_parent", best_pairs$male_parent[1]) data.table::set(results_dt, j, "female_parent", best_pairs$female_parent[1]) data.table::set(results_dt, j, "mendelian_error_pct", base::sprintf("%.2f", min_mismatch_val)) data.table::set(results_dt, j, "markers_tested", base::as.integer(best_markers)) data.table::set(results_dt, j, "status", a_status) + data.table::set(results_dt, j, "tied_candidates", + base::as.integer(base::sum(percent_mismatch == min_mismatch_val, + na.rm = TRUE))) - if (show_ties && num_to_report > 1) { + if (report_all && num_to_report > 1) { tie_row <- base::list(id = prog_id) for (k in base::seq(2, num_to_report)) { tie_row[[base::paste0("male_parent_", k)]] <- best_pairs$male_parent[k] tie_row[[base::paste0("female_parent_", k)]] <- best_pairs$female_parent[k] - tie_row[[base::paste0("mendelian_error_pct_", k)]] <- min_mismatch_val + tie_row[[base::paste0("mendelian_error_pct_", k)]] <- + if (top_n > 1) base::round(percent_mismatch[best_indices[k]], 2) else min_mismatch_val tie_row[[base::paste0("markers_tested_", k)]] <- comparisons[best_indices[k]] } results_list[[j]] <- data.table::as.data.table(tie_row) @@ -395,7 +592,7 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, total_progeny <- base::nrow(final_df) base::cat("\n=== Parentage Assignment Report ===\n") base::cat("\nTotal progeny evaluated:", total_progeny, "\n") - base::cat("Method:", method, " | ", + base::cat("Method:", method_reported, " | ", "Error threshold:", error_threshold, "% | ", "Minimum markers:", min_markers, "\n") diff --git a/R/validate_pedigree.R b/R/validate_pedigree.R index 38a7c12..cf1be87 100644 --- a/R/validate_pedigree.R +++ b/R/validate_pedigree.R @@ -1,8 +1,8 @@ #' Validate Pedigree Trios Using Mendelian Error Analysis #' #' Validates parent-offspring trios against SNP genotype data using Mendelian -#' error rates. Identifies incorrect parentage assignments, suggests -#' best-matching replacements, and outputs a corrected pedigree. Founder trios +#' error rates. Identifies incorrect parentage assignments and outputs a +#' corrected pedigree. Founder trios #' (both parents coded as 0) are preserved unchanged if a founders file is #' supplied. Trios absent from the genotype file are retained as #' no_genotype_data. @@ -47,6 +47,11 @@ #' \item{missing_parents}{Trios with one or both parents coded as 0 (non-founders).} #' \item{full_results}{Complete data.table with all trios and all output columns.} #' \item{corrected_pedigree}{Pedigree table after applying recommended corrections.} +#' \item{marker_summary}{Per-marker Mendelian mismatch counts over every trio +#' with both parents known and genotyped: `marker`, `trios_tested`, +#' `trios_mismatch` and `mismatch_pct`, sorted worst first. Markers with a +#' much higher rate than the rest (e.g. paralogous or miscalled markers) +#' are candidates for removal before rerunning.} #' \item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} #' } #' @@ -120,6 +125,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, ## silence R CMD check NOTEs id <- male_parent <- female_parent <- status <- trio_mendelian_error_pct <- NULL plot_status <- recommended_correction <- NULL + marker <- trios_tested <- trios_mismatch <- mismatch_pct <- NULL #### Input validation #### if (trio_error_threshold < 0 || trio_error_threshold > 100) @@ -200,25 +206,10 @@ validate_pedigree <- function(pedigree_file, genotypes_file, if (base::nrow(pedigree) == 0) stop("No valid trios remain after filtering for genotype availability.") - #### Find best matching parent via homozygous mismatch #### - find_best_parent <- function(prog_id, exclude_ids = base::character(0)) { - candidates <- base::setdiff(base::rownames(genos_hom_mat), - c(prog_id, exclude_ids)) - if (base::length(candidates) == 0) - return(base::list(id = NA_character_, error_pct = NA_real_)) - prog_hom <- genos_hom_mat[prog_id, ] - errors <- base::sapply(candidates, function(cand_id) { - cand_hom <- genos_hom_mat[cand_id, ] - comparisons <- base::sum(!base::is.na(cand_hom) & !base::is.na(prog_hom)) - if (comparisons == 0) return(NA_real_) - (base::sum(cand_hom != prog_hom, na.rm = TRUE) / comparisons) * 100 - }) - if (base::all(base::is.na(errors))) - return(base::list(id = NA_character_, error_pct = NA_real_)) - best_idx <- base::which.min(errors) - base::list(id = candidates[best_idx], - error_pct = base::round(errors[best_idx], 2)) - } + #### Per-marker mismatch counters (filled in the trio loop) #### + marker_names <- base::colnames(genos_mat) + marker_tested_n <- base::numeric(base::length(marker_names)) + marker_mismatch_n <- base::numeric(base::length(marker_names)) #### Main trio evaluation loop #### results_list <- base::lapply(base::seq_len(base::nrow(pedigree)), function(i) { @@ -231,10 +222,6 @@ validate_pedigree <- function(pedigree_file, genotypes_file, markers_tested <- 0L male_parent_error_pct <- NA_real_ female_parent_error_pct <- NA_real_ - best_male_parent <- NA_character_ - best_male_parent_pct <- NA_real_ - best_female_parent <- NA_character_ - best_female_parent_pct <- NA_real_ if (male_parent_id == "0" && female_parent_id == "0" && prog_id %in% founder_ids) { @@ -244,37 +231,25 @@ validate_pedigree <- function(pedigree_file, genotypes_file, if (male_parent_id == "0" && female_parent_id == "0") { status <- "missing_both_parents" correction_decision <- "none" - best_m <- find_best_parent(prog_id, exclude_ids = character(0)) - best_male_parent <- best_m$id - best_male_parent_pct <- best_m$error_pct - best_f <- find_best_parent(prog_id, exclude_ids = c(best_m$id)) - best_female_parent <- best_f$id - best_female_parent_pct <- best_f$error_pct } else if (male_parent_id == "0" && female_parent_id != "0") { status <- "missing_male_parent" correction_decision <- "none" - best_m <- find_best_parent(prog_id, exclude_ids = c(female_parent_id)) - best_male_parent <- best_m$id - best_male_parent_pct <- best_m$error_pct } else if (male_parent_id != "0" && female_parent_id == "0") { status <- "missing_female_parent" correction_decision <- "none" - best_f <- find_best_parent(prog_id, exclude_ids = c(male_parent_id)) - best_female_parent <- best_f$id - best_female_parent_pct <- best_f$error_pct } else { progeny_vec <- genos_mat[prog_id, ] male_parent_vec <- genos_mat[male_parent_id, ] female_parent_vec <- genos_mat[female_parent_id, ] - mismatches <- base::sum( - .mend_mismatch(male_parent_vec, female_parent_vec, - progeny_vec, ploidy), - na.rm = TRUE - ) - markers_tested <- base::sum( - .mend_testable(male_parent_vec, female_parent_vec, - progeny_vec, ploidy) - ) + mm_vec <- .mend_mismatch(male_parent_vec, female_parent_vec, + progeny_vec, ploidy) + tt_vec <- .mend_testable(male_parent_vec, female_parent_vec, + progeny_vec, ploidy) + mismatches <- base::sum(mm_vec, na.rm = TRUE) + markers_tested <- base::sum(tt_vec) + marker_tested_n <<- marker_tested_n + base::as.numeric(tt_vec) + marker_mismatch_n <<- marker_mismatch_n + + base::as.numeric((mm_vec %in% TRUE) & tt_vec) if (markers_tested == 0) { status <- "no_data" correction_decision <- "none" @@ -314,22 +289,10 @@ validate_pedigree <- function(pedigree_file, genotypes_file, correction_decision <- if (error_pct > trio_error_threshold) "unresolved" else "keep_both" } else if (male_acceptable && !female_acceptable) { correction_decision <- "remove_female_parent" - best_f <- find_best_parent(prog_id, exclude_ids = c(male_parent_id)) - best_female_parent <- best_f$id - best_female_parent_pct <- best_f$error_pct } else if (!male_acceptable && female_acceptable) { correction_decision <- "remove_male_parent" - best_m <- find_best_parent(prog_id, exclude_ids = c(female_parent_id)) - best_male_parent <- best_m$id - best_male_parent_pct <- best_m$error_pct } else { correction_decision <- "remove_both" - best_m <- find_best_parent(prog_id, exclude_ids = character(0)) - best_male_parent <- best_m$id - best_male_parent_pct <- best_m$error_pct - best_f <- find_best_parent(prog_id, exclude_ids = c(best_m$id)) - best_female_parent <- best_f$id - best_female_parent_pct <- best_f$error_pct } if (status == "low_markers") correction_decision <- paste0("low_markers_", correction_decision) @@ -347,11 +310,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, status = status, recommended_correction = correction_decision, male_parent_hom_error_pct = male_parent_error_pct, - female_parent_hom_error_pct = female_parent_error_pct, - best_male_candidate = best_male_parent, - best_male_candidate_error_pct = best_male_parent_pct, - best_female_candidate = best_female_parent, - best_female_candidate_error_pct = best_female_parent_pct + female_parent_hom_error_pct = female_parent_error_pct ) }) @@ -368,11 +327,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, status = "no_genotype_data", recommended_correction = "none", male_parent_hom_error_pct = NA_real_, - female_parent_hom_error_pct = NA_real_, - best_male_candidate = NA_character_, - best_male_candidate_error_pct = NA_real_, - best_female_candidate = NA_character_, - best_female_candidate_error_pct = NA_real_ + female_parent_hom_error_pct = NA_real_ ) final_df <- data.table::rbindlist(list(final_df, no_geno_df)) } @@ -487,6 +442,16 @@ validate_pedigree <- function(pedigree_file, genotypes_file, } } + #### Per-marker mismatch summary #### + marker_summary <- data.table::data.table( + marker = marker_names, + trios_tested = base::as.integer(marker_tested_n), + trios_mismatch = base::as.integer(marker_mismatch_n) + ) + marker_summary[, mismatch_pct := base::ifelse( + trios_tested > 0, base::round(trios_mismatch / trios_tested * 100, 2), NA_real_)] + marker_summary <- marker_summary[base::order(-mismatch_pct, -trios_mismatch, marker)] + #### Compile and return named list #### output_list <- base::list( pass = final_df[status == "pass"], @@ -499,6 +464,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, "missing_female_parent")], full_results = final_df, corrected_pedigree = corrected_pedigree, + marker_summary = marker_summary, plot = p ) return(base::invisible(output_list)) diff --git a/README.md b/README.md index ace0952..e54e890 100644 --- a/README.md +++ b/README.md @@ -41,8 +41,8 @@ library(BIGpopA) | Function | Purpose | |---|---| | `check_ped()` | Detect and correct pedigree errors (duplicates, conflicting trios, missing parents, cycles, inconsistent sex roles) | -| `validate_pedigree()` | Validate parent-offspring trios against SNP genotypes using Mendelian error rates | -| `find_parentage()` | Assign the most likely parent(s) to progeny from candidate parents | +| `validate_pedigree()` | Validate parent-offspring trios against SNP genotypes using Mendelian error rates; also returns per-marker mismatch counts (`marker_summary`) | +| `find_parentage()` | Assign the most likely parent(s) to progeny from candidate parents; `top_n` reports several candidates, `tied_candidates` flags ambiguity, and `method = "fill_pedigree"` fills missing parents of a pedigree | | `allele_freq_poly()` | Compute reference population allele frequencies | | `solve_composition_poly()` | Estimate genome-wide breed/line composition by quadratic programming | | `vcf_to_dosage()` | Convert a VCF (`.vcf` / `.vcf.gz`) to allele dosages for any ploidy | diff --git a/cran-comments.md b/cran-comments.md index 256213c..546f191 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -2,24 +2,21 @@ 0 errors | 0 warnings | 0 notes -## Changes in this version (2.1.0) +## Changes in this version (2.2.0) This is a feature update. In this version I have: -* Added `vcf_to_dosage()` to convert VCF files (any ploidy) and `ped_to_dosage()` - to convert 'PLINK' .ped files (diploid) into allele dosages. -* Allowed `find_parentage()`, `validate_pedigree()`, `allele_freq_poly()` and - `solve_composition_poly()` to accept genotypes as text files, VCF files, - 'PLINK' .ped files, data frames or matrices. Existing inputs work unchanged. -* Removed the `ped`, `groups`, `mia`, `sire` and `dam` arguments from - `solve_composition_poly()`. They depended on internal helpers that were not - included in the package and always returned an error, so no working code - is affected. -* Made the `id`, `male_parent`, `female_parent` and `sex` column names - case-insensitive in `check_ped()`, `validate_pedigree()`, `find_parentage()` - and genotype tables. -* Fixed `solve_composition_poly()` failing when only one individual is supplied. -* Examples and tests now remove the temporary files they create. +* Added `top_n` to `find_parentage()` to report the best `top_n` candidates per + progeny, and a `tied_candidates` column that flags ambiguous assignments. + The defaults keep the previous output. +* Added `method = "fill_pedigree"` to `find_parentage()` to fill in missing + parents of a pedigree, with a new `founders_file` argument. +* Added a `marker_summary` element to the result of `validate_pedigree()` with + per-marker Mendelian mismatch counts. +* Removed the suggested replacement parent columns (`best_male_candidate`, `best_male_candidate_error_pct`, `best_female_candidate`, `best_female_candidate_error_pct`) from the `validate_pedigree()` output; `find_parentage(method = "fill_pedigree")` now covers this. Familia does not use these columns. +* Added tests, documentation and vignette sections for the new features. + +Apart from the removed columns, existing calls and outputs are unchanged. 'vcfR' is in Suggests and is only used by `vcf_to_dosage()`; its examples and tests are skipped when it is not installed. @@ -32,5 +29,5 @@ CRAN and dev versions of this package. * We saw 0 new problems * We failed to check 0 packages -Familia only calls `solve_composition_poly(Y, X, ploidy)` and does not use the -removed arguments. +The new arguments and result elements are additive, so Familia's existing +calls are unaffected. diff --git a/man/find_parentage.Rd b/man/find_parentage.Rd index c61ffef..2302d35 100644 --- a/man/find_parentage.Rd +++ b/man/find_parentage.Rd @@ -16,7 +16,9 @@ find_parentage( exclude_self_match = TRUE, verbose = TRUE, plot_results = TRUE, - ploidy = 2 + ploidy = 2, + top_n = 1, + founders_file = NULL ) } \arguments{ @@ -32,10 +34,13 @@ data.table with an 'id' column and an optional 'sex' column ('M', 'F', or 'A'). If absent, all parents are treated as ambiguous.} \item{progeny_file}{Path to a TSV/CSV/TXT file, OR a data.frame / -data.table with an 'id' column.} +data.table with an 'id' column. For \code{method = "fill_pedigree"} it must +be a pedigree with \code{id}, \code{male_parent} and \code{female_parent} columns (e.g. +the \code{corrected_pedigree} from \code{\link[=validate_pedigree]{validate_pedigree()}}); see Details.} \item{method}{Character. One of "best_male_parent", "best_female_parent", -"best_match", or "best_pair" (default).} +"best_match", "best_pair" (default), or "fill_pedigree" (requires a +pedigree as \code{progeny_file}; see Details).} \item{min_markers}{Integer. Minimum markers required; fewer flags low_markers (default: 10).} @@ -66,6 +71,39 @@ test (assumes autopolyploid inheritance; conservative for allopolyploids). Odd ploidy (e.g. triploid), where balanced gametes are undefined, falls back to a model-free opposite-homozygote exclusion evaluated on homozygous-informative markers only (reduced power). Default is 2.} + +\item{founders_file}{Character, optional. Path to a one-column file of +founder IDs. Only used with \code{method = "fill_pedigree"}: founders with +both parents unknown are skipped instead of searched. Default is NULL.} + +\item{top_n}{Integer >= 1. Number of best candidates to report per progeny. +With the default of 1, only the best assignment is reported (plus tied +best pairs when \code{show_ties = TRUE} and \code{method = "best_pair"}). With +\code{top_n > 1}, candidates are ranked by error (ties broken by the number of +markers tested) and the best \code{top_n} are reported: rank 1 in the base +columns and ranks 2 to \code{top_n} in suffix columns (\code{male_parent_2}, +\code{female_parent_2}, \code{mendelian_error_pct_2}, \code{markers_tested_2} for +"best_pair"; \code{best_match_2}, \code{mendelian_error_pct_2}, \code{markers_tested_2} +for the other methods). Ties are then ranked like any other candidate, so +\code{show_ties} has no effect. Progeny with fewer than \code{top_n} candidates +have \code{NA} in the remaining ranks. Status refers to rank 1 only.} +} +\details{ +\strong{Method "fill_pedigree" ("fill in the blanks").} The +\code{progeny_file} must have \code{male_parent} and \code{female_parent} columns. A parent +coded \code{0}, \code{NA} or empty is treated as unknown and each row is handled on +its own: a known male searches only females, a known female searches only +males, and a row with both unknown gets the full pair search (unless listed +in \code{founders_file}). Rows with both parents known are skipped. The search is +the "best_pair" Mendelian trio test with the known parent held fixed, and +candidates come from \code{parents_file} (sex M/F/A as before); the known parent +only needs genotypes. Two columns are added to the results: \code{search_mode} +("pair", "known_male" or "known_female") and \code{known_parent_error_pct} +(homozygous mismatch between the known parent and the progeny; a high value +suggests the known parent is wrong). A known parent absent from the genotype +file gives \code{status = "known_parent_no_genotype"}. Without +\code{method = "fill_pedigree"}, any \code{male_parent} / \code{female_parent} columns in +\code{progeny_file} are ignored. } \value{ A named list (returned invisibly) with elements: @@ -73,7 +111,10 @@ A named list (returned invisibly) with elements: \item{pass}{Progeny with a confident parentage assignment.} \item{high_error}{Progeny whose best assignment exceeds the error threshold.} \item{low_markers}{Progeny with insufficient markers for a valid assignment.} -\item{full_results}{Complete data.table with all progeny and all output columns.} +\item{full_results}{Complete data.table with all progeny and all output columns. +\code{tied_candidates} gives the number of candidates (pairs for "best_pair") +that share the lowest error for each progeny; values above 1 mean the +best assignment is ambiguous.} \item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} } } diff --git a/man/validate_pedigree.Rd b/man/validate_pedigree.Rd index 85655f1..35d5bcf 100644 --- a/man/validate_pedigree.Rd +++ b/man/validate_pedigree.Rd @@ -66,13 +66,18 @@ An invisible named list with the following elements: \item{missing_parents}{Trios with one or both parents coded as 0 (non-founders).} \item{full_results}{Complete data.table with all trios and all output columns.} \item{corrected_pedigree}{Pedigree table after applying recommended corrections.} +\item{marker_summary}{Per-marker Mendelian mismatch counts over every trio +with both parents known and genotyped: \code{marker}, \code{trios_tested}, +\code{trios_mismatch} and \code{mismatch_pct}, sorted worst first. Markers with a +much higher rate than the rest (e.g. paralogous or miscalled markers) +are candidates for removal before rerunning.} \item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} } } \description{ Validates parent-offspring trios against SNP genotype data using Mendelian -error rates. Identifies incorrect parentage assignments, suggests -best-matching replacements, and outputs a corrected pedigree. Founder trios +error rates. Identifies incorrect parentage assignments and outputs a +corrected pedigree. Founder trios (both parents coded as 0) are preserved unchanged if a founders file is supplied. Trios absent from the genotype file are retained as no_genotype_data. diff --git a/tests/testthat/test-find_parentage.R b/tests/testthat/test-find_parentage.R index b2c384d..9a217f1 100644 --- a/tests/testthat/test-find_parentage.R +++ b/tests/testthat/test-find_parentage.R @@ -760,3 +760,151 @@ test_that("invalid input type raises a descriptive error for find_parentage", { regexp = "Error reading input files" ) }) + +# ============================================================================== +# 8. top_n +# ============================================================================== + +test_that("top_n = 3 reports ranked best_pair candidates in suffix columns", { + f <- make_files(base_genos, base_parents, child1_progeny) + out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", top_n = 3, + verbose = FALSE, plot_results = FALSE) + res <- out$full_results + expect_equal(res$male_parent[1], "S1") + expect_equal(res$female_parent[1], "D1") + for (k in 2:3) { + expect_true(all(paste0(c("male_parent_", "female_parent_", + "mendelian_error_pct_", "markers_tested_"), k) %in% names(res))) + expect_false(is.na(res[[paste0("male_parent_", k)]][1])) + } + expect_false("male_parent_4" %in% names(res)) + expect_true(res[["mendelian_error_pct_2"]][1] >= as.numeric(res$mendelian_error_pct[1])) + expect_true(res[["mendelian_error_pct_3"]][1] >= res[["mendelian_error_pct_2"]][1]) +}) + +test_that("top_n = 1 keeps the previous output (no suffix columns without ties)", { + f <- make_files(base_genos, base_parents, child1_progeny) + out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", top_n = 1, + verbose = FALSE, plot_results = FALSE) + expect_false(any(grepl("_\\d+$", names(out$full_results)))) +}) + +test_that("top_n larger than the number of candidate pairs reports all pairs", { + f <- make_files(base_genos, base_parents, child1_progeny) + out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", top_n = 10, + verbose = FALSE, plot_results = FALSE) + expect_true("male_parent_4" %in% names(out$full_results)) + expect_false("male_parent_5" %in% names(out$full_results)) +}) + +test_that("top_n > 1 ranks tied pairs and does not warn when show_ties = FALSE", { + f <- make_files(tied_genos, tied_parents, tied_progeny) + expect_warning( + out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", top_n = 2, + show_ties = FALSE, verbose = FALSE, plot_results = FALSE), + regexp = NA + ) + expect_true("male_parent_2" %in% names(out$full_results)) + expect_false("male_parent_3" %in% names(out$full_results)) +}) + +test_that("top_n works for single-parent methods", { + f <- make_files(base_genos, base_parents, child1_progeny) + out <- find_parentage(f$g, f$p, f$pr, method = "best_male_parent", top_n = 2, + verbose = FALSE, plot_results = FALSE) + res <- out$full_results + expect_equal(res$best_match[1], "S1") + expect_equal(res$best_match_2[1], "S2") + expect_true(res$mendelian_error_pct_2[1] > res$mendelian_error_pct[1]) + expect_true("markers_tested_2" %in% names(res)) +}) + +test_that("top_n keeps progeny order for single-parent methods", { + f <- make_files(base_genos, base_parents, base_progeny) + out <- find_parentage(f$g, f$p, f$pr, method = "best_match", top_n = 2, + verbose = FALSE, plot_results = FALSE) + expect_equal(out$full_results$id, c("child1", "child2")) + expect_equal(out$full_results$best_match, c("S1", "S2")) +}) + +test_that("invalid top_n values are rejected", { + f <- make_files(base_genos, base_parents, child1_progeny) + expect_error(find_parentage(f$g, f$p, f$pr, top_n = 0, verbose = FALSE, plot_results = FALSE), "top_n") + expect_error(find_parentage(f$g, f$p, f$pr, top_n = 1.5, verbose = FALSE, plot_results = FALSE), "top_n") + expect_error(find_parentage(f$g, f$p, f$pr, top_n = "a", verbose = FALSE, plot_results = FALSE), "top_n") +}) + +test_that("tied_candidates counts candidates sharing the lowest error", { + f <- make_files(tied_genos, tied_parents, tied_progeny) + out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", top_n = 2, + verbose = FALSE, plot_results = FALSE) + expect_equal(out$full_results$tied_candidates, 4L) # 2 males x 2 females + f2 <- make_files(base_genos, base_parents, child1_progeny) + out2 <- find_parentage(f2$g, f2$p, f2$pr, method = "best_pair", + verbose = FALSE, plot_results = FALSE) + expect_equal(out2$full_results$tied_candidates, 1L) + out3 <- find_parentage(f2$g, f2$p, f2$pr, method = "best_match", + verbose = FALSE, plot_results = FALSE) + expect_equal(out3$full_results$tied_candidates, 2L) # S1 and D1 both 0% +}) + +# ============================================================================== +# 9. fill_pedigree +# ============================================================================== + +fill_ped <- data.table::data.table( + id = c("child1", "child2", "child2b", "S1", "child_done"), + male_parent = c("S1", "0", "0", "0", "S2"), + female_parent = c("0", "D2", "0", "0", "D2") +) +fill_genos <- rbind(base_genos, + data.table::data.table(id = c("child_done", "child2b"), + M1 = c(2L, 0L), M2 = c(2L, 0L), M3 = c(2L, 0L), + M4 = c(2L, 0L), M5 = c(2L, 0L))) + +test_that("fill_pedigree requires male_parent and female_parent columns", { + f <- make_files(base_genos, base_parents, child1_progeny) + expect_error(find_parentage(f$g, f$p, f$pr, method = "fill_pedigree", + verbose = FALSE, plot_results = FALSE), "fill_pedigree") +}) + +test_that("fill_pedigree fills the missing parent and skips complete trios and founders", { + f <- make_files(fill_genos, base_parents, fill_ped) + fo <- file.path(tempdir(), "founders_fill.txt") + writeLines("S1", fo) + out <- find_parentage(f$g, f$p, f$pr, method = "fill_pedigree", founders_file = fo, + verbose = FALSE, plot_results = FALSE) + res <- as.data.frame(out$full_results) + expect_false("child_done" %in% res$id) # complete trio + expect_false("S1" %in% res$id) # founder + expect_equal(res$female_parent[res$id == "child1"], "D1") + expect_equal(res$male_parent[res$id == "child1"], "S1") + expect_equal(res$search_mode[res$id == "child1"], "known_male") + expect_equal(res$male_parent[res$id == "child2"], "S2") + expect_equal(res$female_parent[res$id == "child2"], "D2") + expect_equal(res$search_mode[res$id == "child2"], "known_female") + expect_equal(res$search_mode[res$id == "child2b"], "pair") + expect_true(all(c("known_parent_error_pct", "tied_candidates") %in% names(res))) +}) + +test_that("fill_pedigree restricts candidates to the opposite sex and reports top_n", { + f <- make_files(base_genos, base_parents, data.table::data.table( + id = "child1", male_parent = "S1", female_parent = "0")) + out <- find_parentage(f$g, f$p, f$pr, method = "fill_pedigree", top_n = 2, + verbose = FALSE, plot_results = FALSE) + res <- out$full_results + expect_equal(res$male_parent[1], "S1") + expect_equal(res$female_parent[1], "D1") # all-0 child: only D1 fits S1 (all 0) + expect_equal(res$male_parent_2[1], "S1") # known parent stays fixed + expect_equal(res$female_parent_2[1], "D2") +}) + +test_that("fill_pedigree flags a known parent without genotypes", { + f <- make_files(base_genos, base_parents, data.table::data.table( + id = "child1", male_parent = "ghost", female_parent = "0")) + expect_warning( + out <- find_parentage(f$g, f$p, f$pr, method = "fill_pedigree", + verbose = FALSE, plot_results = FALSE), + "ghost") + expect_equal(out$full_results$status[1], "known_parent_no_genotype") +}) diff --git a/tests/testthat/test-validate_pedigree.R b/tests/testthat/test-validate_pedigree.R index 520bd62..38c6288 100644 --- a/tests/testthat/test-validate_pedigree.R +++ b/tests/testthat/test-validate_pedigree.R @@ -181,9 +181,7 @@ test_that("full_results has all expected lowercase columns", { "id", "orig_male_parent", "orig_female_parent", "trio_mendelian_error_pct", "trio_markers_tested", "status", "recommended_correction", - "male_parent_hom_error_pct", "female_parent_hom_error_pct", - "best_male_candidate", "best_male_candidate_error_pct", - "best_female_candidate", "best_female_candidate_error_pct" + "male_parent_hom_error_pct", "female_parent_hom_error_pct" ) expect_true(all(expected_cols %in% names(out$full_results))) }) @@ -366,8 +364,6 @@ test_that("missing_male_parent status and recommendation are correct", { r <- out$full_results[id == "IND_E"] expect_equal(r$status, "missing_male_parent") expect_equal(r$recommended_correction, "none") - expect_false(is.na(r$best_male_candidate)) - expect_true(is.na(r$best_female_candidate)) }) test_that("missing_female_parent status and recommendation are correct", { @@ -379,8 +375,6 @@ test_that("missing_female_parent status and recommendation are correct", { r <- out$full_results[id == "IND_E"] expect_equal(r$status, "missing_female_parent") expect_equal(r$recommended_correction, "none") - expect_true(is.na(r$best_male_candidate)) - expect_false(is.na(r$best_female_candidate)) }) test_that("missing_both_parents status and recommendations are correct", { @@ -392,28 +386,6 @@ test_that("missing_both_parents status and recommendations are correct", { r <- out$full_results[id == "IND_E"] expect_equal(r$status, "missing_both_parents") expect_equal(r$recommended_correction, "none") - expect_false(is.na(r$best_male_candidate)) - expect_false(is.na(r$best_female_candidate)) -}) - -test_that("best_male_candidate for missing_male_parent excludes the known female parent", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_E", male_parent = "0", - female_parent = "IND_B")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_E"] - expect_false(r$best_male_candidate == "IND_B") -}) - -test_that("best_female_candidate for missing_female_parent excludes the known male parent", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_E", male_parent = "IND_A", - female_parent = "0")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_E"] - expect_false(r$best_female_candidate == "IND_A") }) test_that("missing_parents list element contains only missing_* rows", { @@ -443,8 +415,6 @@ test_that("founders status is assigned when ID is in founders list with 0 0 pare r <- out$full_results[id == "IND_A"] expect_equal(r$status, "founders") expect_equal(r$recommended_correction, "none") - expect_true(is.na(r$best_male_candidate)) - expect_true(is.na(r$best_female_candidate)) }) test_that("founders list element contains only founders rows", { @@ -524,8 +494,6 @@ test_that("no_genotype_data rows have NA/0 for all analysis columns", { r <- out$full_results[id == "GHOST"] expect_true(is.na(r$trio_mendelian_error_pct)) expect_equal(r$trio_markers_tested, 0L) - expect_true(is.na(r$best_male_candidate)) - expect_true(is.na(r$best_female_candidate)) }) test_that("no_genotype_data flagged when a declared parent is absent from genotype file", { diff --git a/tests/testthat/test-validate_pedigree_marker_summary.R b/tests/testthat/test-validate_pedigree_marker_summary.R new file mode 100644 index 0000000..b9f765f --- /dev/null +++ b/tests/testthat/test-validate_pedigree_marker_summary.R @@ -0,0 +1,113 @@ +# tests/testthat/test-validate_pedigree_marker_summary.R +# Independent check of validate_pedigree()$marker_summary against a brute-force +# per-marker oracle (gamete-range rule), on simulated trios with known errors. + +library(testthat) +library(data.table) + +ms_sim <- function(ploidy, nm = 200, seed = 7) { + set.seed(seed + ploidy) + p <- runif(nm, 0.2, 0.8) + draw <- function() rbinom(nm, ploidy, p) + gam <- function(g) rhyper(length(g), g, ploidy - g, ploidy / 2) + G <- list() + for (id in c("M1", "M2", "M3", "F1", "F2", "F3")) G[[id]] <- draw() + cross <- function(m, f) gam(G[[m]]) + gam(G[[f]]) + G[["C1"]] <- cross("M1", "F1") + c2 <- cross("M2", "F2"); c2[sample(nm, 6)] <- NA # missing calls + G[["C2"]] <- c2 + c3 <- cross("M3", "F3"); idx <- sample(nm, 10); c3[idx] <- draw()[idx] + G[["C3"]] <- c3 # genotyping errors + G[["C4"]] <- cross("M1", "F2") # clean + G[["C5"]] <- cross("M2", "F3") # parent F3 genotyped + mat <- do.call(rbind, G) + colnames(mat) <- paste0("m", seq_len(nm)) + geno <- data.table(id = rownames(mat), as.data.table(mat)) + ped <- data.table( + id = c("M1", "M2", "M3", "F1", "F2", "F3", "C1", "C2", "C3", "C4", "C5", "C6", "C7"), + male_parent = c(NA, NA, NA, NA, NA, NA, "M1", "M2", "M3", "M1", "M2", "M1", "M1"), + female_parent = c(NA, NA, NA, NA, NA, NA, "F1", "F2", "F3", "F2", "F3", NA, "FX")) + # C6: one parent unknown; C7: female parent "FX" has no genotype -> neither is a full trio + list(geno = geno, ped = ped, GM = mat, ploidy = ploidy) +} + +ms_oracle <- function(pop) { + ped <- as.data.frame(pop$ped); GM <- pop$GM; h <- pop$ploidy / 2 + keep <- !is.na(ped$male_parent) & !is.na(ped$female_parent) & + ped$id %in% rownames(GM) & ped$male_parent %in% rownames(GM) & + ped$female_parent %in% rownames(GM) + ped <- ped[keep, , drop = FALSE] + tested <- mism <- numeric(ncol(GM)) + for (i in seq_len(nrow(ped))) { + m <- GM[ped$male_parent[i], ]; f <- GM[ped$female_parent[i], ]; o <- GM[ped$id[i], ] + ok <- !is.na(m) & !is.na(f) & !is.na(o) + lo <- pmax(0, m - h) + pmax(0, f - h) + hi <- pmin(m, h) + pmin(f, h) + bad <- ok & (o < lo | o > hi) + tested <- tested + ok + mism <- mism + bad + } + data.frame(marker = colnames(GM), trios_tested = tested, trios_mismatch = mism, + stringsAsFactors = FALSE) +} + +run_ms <- function(pop) + suppressWarnings(validate_pedigree(pop$ped, pop$geno, min_markers = 10, + verbose = FALSE, plot_results = FALSE, + ploidy = pop$ploidy)) + +for (pl in c(2, 4)) { + test_that(paste0("marker_summary counts match oracle (ploidy ", pl, ")"), { + pop <- ms_sim(pl) + res <- run_ms(pop) + ms <- as.data.frame(res$marker_summary) + ex <- ms_oracle(pop) + + expect_true(all(c("marker", "trios_tested", "trios_mismatch", "mismatch_pct") %in% names(ms))) + expect_equal(nrow(ms), ncol(pop$GM)) + expect_setequal(ms$marker, ex$marker) + + m <- merge(ms, ex, by = "marker", suffixes = c("", ".ex")) + expect_equal(as.numeric(m$trios_tested), as.numeric(m$trios_tested.ex)) + expect_equal(as.numeric(m$trios_mismatch), as.numeric(m$trios_mismatch.ex)) + + exp_pct <- ifelse(m$trios_tested.ex > 0, + round(m$trios_mismatch.ex / m$trios_tested.ex * 100, 2), NA) + expect_equal(m$mismatch_pct, exp_pct) + + # mismatches never exceed tested; at most the number of full trios + expect_true(all(ms$trios_mismatch <= ms$trios_tested)) + expect_true(max(ms$trios_tested) <= 5) # C1-C5 only + }) + + test_that(paste0("marker_summary is sorted by mismatch_pct desc (ploidy ", pl, ")"), { + ms <- as.data.frame(run_ms(ms_sim(pl))$marker_summary) + pct <- ms$mismatch_pct + pct[is.na(pct)] <- -Inf + expect_false(is.unsorted(rev(pct))) + # ties on pct broken by mismatch count desc + same <- which(diff(pct) == 0) + expect_true(all(diff(ms$trios_mismatch)[same] <= 0)) + }) +} + +test_that("a clean population gives zero mismatches everywhere", { + pop <- ms_sim(2) + pop$ped <- pop$ped[is.na(male_parent) | id %in% c("C1", "C4")] + ms <- as.data.frame(run_ms(pop)$marker_summary) + expect_equal(sum(ms$trios_mismatch), 0) + expect_equal(sum(ms_oracle(pop)$trios_mismatch), 0) +}) + +test_that("a known injected error is attributed to the right marker", { + pop <- ms_sim(2) + g <- data.table::copy(pop$geno) + g[id == "M1", m1 := 0L]; g[id == "F1", m1 := 0L]; g[id == "C1", m1 := 2L] + pop$geno <- g + pop$GM["M1", "m1"] <- 0; pop$GM["F1", "m1"] <- 0; pop$GM["C1", "m1"] <- 2 + pop$ped <- pop$ped[is.na(male_parent) | id == "C1"] + ms <- as.data.frame(run_ms(pop)$marker_summary) + expect_equal(ms$trios_mismatch[ms$marker == "m1"], 1) + expect_equal(ms$trios_tested[ms$marker == "m1"], 1) + expect_equal(ms$mismatch_pct[ms$marker == "m1"], 100) +}) diff --git a/vignettes/BIGpopA.Rmd b/vignettes/BIGpopA.Rmd index 4f81ffc..cbe3785 100644 --- a/vignettes/BIGpopA.Rmd +++ b/vignettes/BIGpopA.Rmd @@ -35,9 +35,10 @@ functions in sequence: * check_ped() — detect and correct structural pedigree errors (duplicates, conflicting trios, inconsistent sex roles, missing parents). * validate_pedigree() — test parent–offspring trios for Mendelian - errors against SNP genotypes and flag incorrect assignments. + errors against SNP genotypes, flag incorrect assignments and summarize + mismatches per marker. * find_parentage() — assign the most likely parent(s) to progeny from a - pool of candidate parents. + pool of candidate parents, or fill in the missing parents of a pedigree. It also estimates genome-wide breed/line composition for diploid and polyploid populations with allele_freq_poly() and solve_composition_poly(), shown at the @@ -138,8 +139,8 @@ clean_ped_results$missing_parents. validate_pedigree() compares each trio against the genotype matrix and computes a Mendelian error rate. Trios are classified as pass, fail, low_markers, -no_genotype_data, founders, or missing_parents, and a best-matching replacement -parent is suggested where relevant. With plot_results = TRUE the function draws +no_genotype_data, founders, or missing_parents. To fill in missing or removed parents, +run find_parentage() with method = "fill_pedigree". With plot_results = TRUE the function draws a histogram of trio error percentages. ```{r} @@ -158,6 +159,17 @@ genotyped; and C_missing has no_genotype_data. Key tuning arguments are trio_error_threshold (default 5.0), min_markers (default 10), and single_parent_error_threshold (default 2.0). +The result also contains marker_summary, a per-marker table with the number of +trios tested, the number with a Mendelian mismatch, and the mismatch +percentage, sorted worst first. It covers every trio with both parents known and +genotyped, and helps find markers (for example paralogous or miscalled SNPs) +that drive many failures. + +```{r} +#marker_summary +head(ped_validate_results$marker_summary) +``` + # Parentage assignment When parentage is unknown, find_parentage() searches a pool of candidate @@ -182,6 +194,36 @@ its recorded parents, is here correctly recovered to its true pair (M1/F1) at error_threshold, min_markers, show_ties, allow_parent_selfing, and exclude_self_match. +The tied_candidates column gives the number of candidates (pairs for best_pair) +that share each progeny's lowest error; values above 1 mean the best assignment +is ambiguous. Set top_n to report the best top_n candidates per progeny; rank 1 +stays in the base columns and ranks 2 to top_n are added as suffix columns +(for example male_parent_2, female_parent_2, mendelian_error_pct_2). + +```{r} +#top_n +find_parentage(genotypes, parents, progeny, method = "best_pair", top_n = 3, + verbose = FALSE, plot_results = FALSE)$full_results +``` + +## Filling in missing parents of a pedigree + +With method = "fill_pedigree", progeny is a pedigree with id, male_parent and +female_parent columns (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 file +(with sex M, F or A) is still required. The founders_file argument skips +founders. The output adds search_mode and known_parent_error_pct. + +```{r} +#fill_pedigree +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 +``` # Breed and line composition Beyond pedigree QC, BIGpopA estimates genome-wide breed or line composition for From 95f6290e666498880326b95671096e7ed786b00d Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Thu, 8 Oct 2026 16:54:01 -0400 Subject: [PATCH 4/8] fix validate_pedigree low_markers correction, add assign_composition_line and marker_trio_table --- NAMESPACE | 1 + NEWS.md | 3 + R/assign_composition_line.R | 52 ++++++++++++++++ R/utils.R | 48 +++++++++++++++ R/validate_pedigree.R | 59 +++++++++++++------ README.md | 1 + cran-comments.md | 16 +++++ man/assign_composition_line.Rd | 41 +++++++++++++ man/find_parentage.Rd | 42 ++++++------- man/validate_pedigree.Rd | 21 ++++++- tests/testthat/test-assign_composition_line.R | 42 +++++++++++++ tests/testthat/test-validate_pedigree.R | 32 +++++++++- .../test-validate_pedigree_marker_trio.R | 55 +++++++++++++++++ vignettes/BIGpopA.Rmd | 20 +++---- 14 files changed, 378 insertions(+), 55 deletions(-) create mode 100644 R/assign_composition_line.R create mode 100644 man/assign_composition_line.Rd create mode 100644 tests/testthat/test-assign_composition_line.R create mode 100644 tests/testthat/test-validate_pedigree_marker_trio.R diff --git a/NAMESPACE b/NAMESPACE index 234c7c0..de22791 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,7 @@ # Generated by roxygen2: do not edit by hand export(allele_freq_poly) +export(assign_composition_line) export(check_ped) export(find_parentage) export(ped_to_dosage) diff --git a/NEWS.md b/NEWS.md index 86483ce..66f583c 100644 --- a/NEWS.md +++ b/NEWS.md @@ -4,6 +4,9 @@ - `find_parentage()`: new `method = "fill_pedigree"` fills in missing parents of a pedigree. `progeny_file` must have `id`, `male_parent` and `female_parent` columns (e.g. `validate_pedigree()$corrected_pedigree`); parents coded `0`/`NA`/empty are unknown. A known male searches only females, a known female only males, and rows with both unknown get the full pair search (skip founders with the new `founders_file`). Rows with both parents known are skipped. Results gain `search_mode` and `known_parent_error_pct`, and a known parent without genotypes gets `status = "known_parent_no_genotype"`. - `find_parentage()`: new `top_n` argument reports the best `top_n` candidates per progeny, ranked by Mendelian error (ties broken by markers tested). 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`, `markers_tested_2` for `best_pair`; `best_match_2` and the same error/marker columns for the other methods). The default `top_n = 1` keeps the previous output, including `show_ties`. - `find_parentage()`: new `tied_candidates` column gives the number of candidates (pairs for `best_pair`) that share each progeny's lowest error. Values above 1 mean the best assignment is ambiguous, including for the single-parent methods, which previously did not report ties. +- `validate_pedigree()`: fixed `corrected_pedigree` for `low_markers_remove_female_parent`, which also set the male parent to `0`. Low-marker trios now remove only the identified parent(s). +- New `assign_composition_line()` turns the output of `solve_composition_poly()` into the composition table used by Familia: percentages per reference population and a `Predicted line`, labelled `Undetermined` when the highest proportion is below `assign_threshold` (default 50). `solve_composition_poly()` is unchanged. +- `validate_pedigree()`: new `marker_trio_table = FALSE` argument. When TRUE the result includes `marker_trio_table`, a marker-by-trio table with one column per trio (`id_male_parent_female_parent`) giving each marker's status: `match`, `mismatch_male`, `mismatch_female`, `mismatch_both`, `mismatch_combination` (each parent compatible alone but not the pair), `untestable` (odd ploidy) or `missing_*` named after the individuals without a call (e.g. `missing_male`, `missing_male_female`, `missing_progeny`). # BIGpopA 2.1.1 - `validate_pedigree()`: a trio that fails the Mendelian error threshold while neither parent fails the single-parent homozygous check now gets `recommended_correction = "unresolved"` instead of `"keep_both"`. This happens, for example, with a recorded self whose offspring carries alleles the parent lacks: the pedigree is wrong, but the faulty parent cannot be identified. These trios are left unchanged in `corrected_pedigree` for manual review. diff --git a/R/assign_composition_line.R b/R/assign_composition_line.R new file mode 100644 index 0000000..9b4877a --- /dev/null +++ b/R/assign_composition_line.R @@ -0,0 +1,52 @@ +#' Assign Each Sample to a Reference Line +#' +#' Builds the composition table shown in the Familia BreedTools tab from the +#' output of [solve_composition_poly()]: one row per sample with the estimated +#' percentage of each reference population and the predicted line. A sample is +#' assigned to the population with its highest proportion only if that +#' proportion is at least `assign_threshold`; otherwise it is labelled +#' `"Undetermined"`. +#' +#' @param composition Matrix or data.frame returned by +#' [solve_composition_poly()]: samples in rows, reference populations in +#' columns (proportions summing to 1). An `R2` column, if present, is dropped. +#' @param assign_threshold Numeric between 0 and 100. Minimum percentage the +#' highest population must reach for the sample to be assigned to it +#' (default: 50). Use 0 to always assign the highest population. +#' +#' @return A data.frame with an `ID` column, one column per reference +#' population holding the estimated percentage (0 to 100, rounded to one +#' decimal) and a `Predicted line` column. Ties between populations go to the +#' first column. +#' +#' @examples +#' comp <- rbind( +#' Sample01 = c(Duroc = 0.998, Hampshire = 0.000, Landrace = 0.002, R2 = 0.991), +#' Sample02 = c(Duroc = 0.529, Hampshire = 0.471, Landrace = 0.000, R2 = 0.973), +#' Sample03 = c(Duroc = 0.340, Hampshire = 0.330, Landrace = 0.330, R2 = 0.952) +#' ) +#' assign_composition_line(comp) +#' assign_composition_line(comp, assign_threshold = 0) +#' +#' @export +assign_composition_line <- function(composition, assign_threshold = 50) { + if (!is.numeric(assign_threshold) || length(assign_threshold) != 1 || + is.na(assign_threshold) || assign_threshold < 0 || assign_threshold > 100) + stop("assign_threshold must be a single number between 0 and 100.") + + comp <- as.data.frame(composition, check.names = FALSE) + comp <- comp[, !colnames(comp) %in% "R2", drop = FALSE] + if (ncol(comp) == 0) stop("composition has no population columns.") + comp[] <- lapply(comp, as.numeric) + mat <- as.matrix(comp) + + top_value <- apply(mat, 1, max, na.rm = TRUE) + top_line <- colnames(mat)[max.col(mat, ties.method = "first")] + predicted <- ifelse(top_value * 100 >= assign_threshold, top_line, "Undetermined") + + out <- data.frame(ID = rownames(mat), round(mat * 100, 1), + check.names = FALSE, stringsAsFactors = FALSE, + row.names = NULL) + out[["Predicted line"]] <- predicted + out +} diff --git a/R/utils.R b/R/utils.R index 548241a..6027216 100644 --- a/R/utils.R +++ b/R/utils.R @@ -301,3 +301,51 @@ QPsolve <- function(Y, X) { return(result) } +#' Per-marker status of one trio +#' +#' Classifies every marker of one parent-offspring trio as a match, a mismatch +#' attributed to the male parent, the female parent or both, a mismatch that only +#' appears for the pair of parents (\code{mismatch_combination}), a marker that +#' cannot be tested (\code{untestable}, odd ploidy only) or a missing call (named after every individual without a call, for example \code{missing_male_female}). For +#' even ploidy a parent is incompatible with the offspring when the offspring +#' dose lies outside that parent's gamete range plus any contribution (0 to +#' ploidy / 2 alleles) from the other parent. For odd ploidy a parent is +#' incompatible when it is homozygous for the opposite allele of a homozygous +#' offspring. +#' +#' @param male,female,offspring dosage vectors (0..ploidy), aligned by marker. +#' @param ploidy integer ploidy level. +#' @return Character vector with one status per marker. +#' @noRd +.marker_trio_status <- function(male, female, offspring, ploidy) { + mm <- .mend_mismatch(male, female, offspring, ploidy) + tt <- .mend_testable(male, female, offspring, ploidy) + if (ploidy %% 2 == 0) { + h <- ploidy / 2 + parent_bad <- function(g) { + !(offspring >= base::pmax(0, g - h) & offspring <= base::pmin(g, h) + h) + } + male_bad <- parent_bad(male) + female_bad <- parent_bad(female) + } else { + o_hom <- offspring == 0 | offspring == ploidy + male_bad <- (male == 0 | male == ploidy) & male != offspring & o_hom + female_bad <- (female == 0 | female == ploidy) & female != offspring & o_hom + } + # missing calls are named after every individual without a call, e.g. + # missing_male, missing_male_female, missing_progeny_female + na_p <- base::is.na(offspring); na_m <- base::is.na(male); na_f <- base::is.na(female) + missing_label <- base::paste0("missing", + base::ifelse(na_p, "_progeny", ""), + base::ifelse(na_m, "_male", ""), + base::ifelse(na_f, "_female", "")) + dplyr::case_when( + na_p | na_m | na_f ~ missing_label, + !tt ~ "untestable", + !mm ~ "match", + male_bad & !female_bad ~ "mismatch_male", + female_bad & !male_bad ~ "mismatch_female", + male_bad & female_bad ~ "mismatch_both", + TRUE ~ "mismatch_combination" + ) +} diff --git a/R/validate_pedigree.R b/R/validate_pedigree.R index cf1be87..ea723fb 100644 --- a/R/validate_pedigree.R +++ b/R/validate_pedigree.R @@ -36,6 +36,9 @@ #' Odd ploidy (e.g. triploid), where balanced gametes are undefined, falls #' back to a model-free opposite-homozygote exclusion evaluated on #' homozygous-informative markers only (reduced power). Default is 2. +#' @param marker_trio_table Logical. If TRUE, the returned list includes +#' `marker_trio_table`, a marker-by-trio table of per-marker statuses (see +#' Details). Default is FALSE. #' #' @return An invisible named list with the following elements: #' \describe{ @@ -52,10 +55,24 @@ #' `trios_mismatch` and `mismatch_pct`, sorted worst first. Markers with a #' much higher rate than the rest (e.g. paralogous or miscalled markers) #' are candidates for removal before rerunning.} +#' \item{marker_trio_table}{If `marker_trio_table = TRUE`, a data.table with +#' one row per marker (`marker`) and one column per trio with both parents +#' known and genotyped, named `id_male_parent_female_parent`; otherwise NULL. +#' See Details for the status values.} #' \item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} #' } #' -#' @details The `recommended_correction` column of `full_results` takes these +#' @details The cells of `marker_trio_table` take one of these values: `match`; +#' `mismatch_male` or `mismatch_female` (the offspring is impossible with that +#' parent alone, whatever the other parent contributes); `mismatch_both` (the +#' offspring is impossible with either parent alone); `mismatch_combination` +#' (each parent is compatible alone but the pair is not, for example a 0 x 0 +#' cross with a heterozygous offspring); `untestable` (odd ploidy, marker not +#' informative); `missing_*` (no call, named after every individual without a +#' call: `missing_progeny`, `missing_male`, `missing_female`, `missing_male_female`, +#' `missing_progeny_male`, `missing_progeny_female` or `missing_progeny_male_female`). +#' +#' The `recommended_correction` column of `full_results` takes these #' values: #' \describe{ #' \item{none}{Trio passes, or a parent is missing / no genotype data.} @@ -120,7 +137,8 @@ validate_pedigree <- function(pedigree_file, genotypes_file, single_parent_error_threshold = 2.0, verbose = TRUE, plot_results = TRUE, - ploidy = 2) { + ploidy = 2, + marker_trio_table = FALSE) { ## silence R CMD check NOTEs id <- male_parent <- female_parent <- status <- trio_mendelian_error_pct <- NULL @@ -338,23 +356,12 @@ validate_pedigree <- function(pedigree_file, genotypes_file, prog_id <- final_df$id[i] decision <- final_df$recommended_correction[i] row_idx <- base::which(corrected_pedigree$id == prog_id) - if (decision == "remove_male_parent") { - data.table::set(corrected_pedigree, row_idx, "male_parent", "0") - } else if (decision == "remove_female_parent") { + # low_markers_* decisions apply the same removal as the plain decision + base_decision <- base::sub("^low_markers_", "", decision) + if (base_decision %in% c("remove_male_parent", "remove_both")) + data.table::set(corrected_pedigree, row_idx, "male_parent", "0") + if (base_decision %in% c("remove_female_parent", "remove_both")) data.table::set(corrected_pedigree, row_idx, "female_parent", "0") - } else if (decision %in% c("remove_both", - "low_markers_remove_both", - "low_markers_remove_male_parent", - "low_markers_remove_female_parent")) { - if (grepl("male", decision)) - data.table::set(corrected_pedigree, row_idx, "male_parent", "0") - if (grepl("female", decision)) - data.table::set(corrected_pedigree, row_idx, "female_parent", "0") - if (decision %in% c("low_markers_remove_both", "remove_both")) { - data.table::set(corrected_pedigree, row_idx, "male_parent", "0") - data.table::set(corrected_pedigree, row_idx, "female_parent", "0") - } - } } #### Summary output #### @@ -442,6 +449,21 @@ validate_pedigree <- function(pedigree_file, genotypes_file, } } + #### Optional marker x trio status table #### + marker_trio_out <- NULL + if (marker_trio_table) { + full_trios <- pedigree[male_parent != "0" & female_parent != "0"] + status_list <- base::lapply(base::seq_len(base::nrow(full_trios)), function(i) + .marker_trio_status(genos_mat[full_trios$male_parent[i], ], + genos_mat[full_trios$female_parent[i], ], + genos_mat[full_trios$id[i], ], + ploidy)) + base::names(status_list) <- base::make.unique( + base::paste(full_trios$id, full_trios$male_parent, full_trios$female_parent, sep = "_")) + marker_trio_out <- data.table::as.data.table( + base::c(base::list(marker = marker_names), status_list)) + } + #### Per-marker mismatch summary #### marker_summary <- data.table::data.table( marker = marker_names, @@ -465,6 +487,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, full_results = final_df, corrected_pedigree = corrected_pedigree, marker_summary = marker_summary, + marker_trio_table = marker_trio_out, plot = p ) return(base::invisible(output_list)) diff --git a/README.md b/README.md index e54e890..d84f990 100644 --- a/README.md +++ b/README.md @@ -45,6 +45,7 @@ library(BIGpopA) | `find_parentage()` | Assign the most likely parent(s) to progeny from candidate parents; `top_n` reports several candidates, `tied_candidates` flags ambiguity, and `method = "fill_pedigree"` fills missing parents of a pedigree | | `allele_freq_poly()` | Compute reference population allele frequencies | | `solve_composition_poly()` | Estimate genome-wide breed/line composition by quadratic programming | +| `assign_composition_line()` | Percentages per line and predicted line, with an assignment threshold (default 50%) | | `vcf_to_dosage()` | Convert a VCF (`.vcf` / `.vcf.gz`) to allele dosages for any ploidy | | `ped_to_dosage()` | Convert a PLINK `.ped` (with optional `.map`) to allele dosages | diff --git a/cran-comments.md b/cran-comments.md index 546f191..ca397e9 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -14,10 +14,26 @@ This is a feature update. In this version I have: * Added a `marker_summary` element to the result of `validate_pedigree()` with per-marker Mendelian mismatch counts. * Removed the suggested replacement parent columns (`best_male_candidate`, `best_male_candidate_error_pct`, `best_female_candidate`, `best_female_candidate_error_pct`) from the `validate_pedigree()` output; `find_parentage(method = "fill_pedigree")` now covers this. Familia does not use these columns. +* Added a `marker_trio_table` argument to `validate_pedigree()` that returns a marker-by-trio table of per-marker statuses (default FALSE, so existing output is unchanged). +* Added `assign_composition_line()`, which turns the output of `solve_composition_poly()` into a table of percentages per line with a predicted line and an assignment threshold (default 50%). `solve_composition_poly()` is unchanged. +* Fixed `validate_pedigree()` so that `low_markers_remove_female_parent` removes only the female parent in `corrected_pedigree` (it previously also removed the male parent). * Added tests, documentation and vignette sections for the new features. Apart from the removed columns, existing calls and outputs are unchanged. +## Why this is version 2.2.0 and not 3.0.0 + +The four `best_*_candidate` columns removed from `validate_pedigree()` were +added in 2.1.0 as an experimental suggestion that duplicated, with weaker +statistics, what `find_parentage()` does. They are the only output removed and +no function argument was removed or changed. The only reverse dependency +(Familia) does not use them, as shown by the reverse dependency check below, +and all other calls and outputs are unchanged. Because the break is limited to +four optional output columns with a documented replacement +(`find_parentage(method = "fill_pedigree")`) and has no effect on any package +on CRAN, we judged a minor version increase appropriate. The removal is stated +in `NEWS.md`. + 'vcfR' is in Suggests and is only used by `vcf_to_dosage()`; its examples and tests are skipped when it is not installed. diff --git a/man/assign_composition_line.Rd b/man/assign_composition_line.Rd new file mode 100644 index 0000000..97ceb51 --- /dev/null +++ b/man/assign_composition_line.Rd @@ -0,0 +1,41 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/assign_composition_line.R +\name{assign_composition_line} +\alias{assign_composition_line} +\title{Assign Each Sample to a Reference Line} +\usage{ +assign_composition_line(composition, assign_threshold = 50) +} +\arguments{ +\item{composition}{Matrix or data.frame returned by +\code{\link[=solve_composition_poly]{solve_composition_poly()}}: samples in rows, reference populations in +columns (proportions summing to 1). An \code{R2} column, if present, is dropped.} + +\item{assign_threshold}{Numeric between 0 and 100. Minimum percentage the +highest population must reach for the sample to be assigned to it +(default: 50). Use 0 to always assign the highest population.} +} +\value{ +A data.frame with an \code{ID} column, one column per reference +population holding the estimated percentage (0 to 100, rounded to one +decimal) and a \verb{Predicted line} column. Ties between populations go to the +first column. +} +\description{ +Builds the composition table shown in the Familia BreedTools tab from the +output of \code{\link[=solve_composition_poly]{solve_composition_poly()}}: one row per sample with the estimated +percentage of each reference population and the predicted line. A sample is +assigned to the population with its highest proportion only if that +proportion is at least \code{assign_threshold}; otherwise it is labelled +\code{"Undetermined"}. +} +\examples{ +comp <- rbind( + Sample01 = c(Duroc = 0.998, Hampshire = 0.000, Landrace = 0.002, R2 = 0.991), + Sample02 = c(Duroc = 0.529, Hampshire = 0.471, Landrace = 0.000, R2 = 0.973), + Sample03 = c(Duroc = 0.340, Hampshire = 0.330, Landrace = 0.330, R2 = 0.952) +) +assign_composition_line(comp) +assign_composition_line(comp, assign_threshold = 0) + +} diff --git a/man/find_parentage.Rd b/man/find_parentage.Rd index 2302d35..b8b9b27 100644 --- a/man/find_parentage.Rd +++ b/man/find_parentage.Rd @@ -72,10 +72,6 @@ Odd ploidy (e.g. triploid), where balanced gametes are undefined, falls back to a model-free opposite-homozygote exclusion evaluated on homozygous-informative markers only (reduced power). Default is 2.} -\item{founders_file}{Character, optional. Path to a one-column file of -founder IDs. Only used with \code{method = "fill_pedigree"}: founders with -both parents unknown are skipped instead of searched. Default is NULL.} - \item{top_n}{Integer >= 1. Number of best candidates to report per progeny. With the default of 1, only the best assignment is reported (plus tied best pairs when \code{show_ties = TRUE} and \code{method = "best_pair"}). With @@ -87,23 +83,10 @@ columns and ranks 2 to \code{top_n} in suffix columns (\code{male_parent_2}, for the other methods). Ties are then ranked like any other candidate, so \code{show_ties} has no effect. Progeny with fewer than \code{top_n} candidates have \code{NA} in the remaining ranks. Status refers to rank 1 only.} -} -\details{ -\strong{Method "fill_pedigree" ("fill in the blanks").} The -\code{progeny_file} must have \code{male_parent} and \code{female_parent} columns. A parent -coded \code{0}, \code{NA} or empty is treated as unknown and each row is handled on -its own: a known male searches only females, a known female searches only -males, and a row with both unknown gets the full pair search (unless listed -in \code{founders_file}). Rows with both parents known are skipped. The search is -the "best_pair" Mendelian trio test with the known parent held fixed, and -candidates come from \code{parents_file} (sex M/F/A as before); the known parent -only needs genotypes. Two columns are added to the results: \code{search_mode} -("pair", "known_male" or "known_female") and \code{known_parent_error_pct} -(homozygous mismatch between the known parent and the progeny; a high value -suggests the known parent is wrong). A known parent absent from the genotype -file gives \code{status = "known_parent_no_genotype"}. Without -\code{method = "fill_pedigree"}, any \code{male_parent} / \code{female_parent} columns in -\code{progeny_file} are ignored. + +\item{founders_file}{Character, optional. Path to a one-column file of +founder IDs. Only used when \code{progeny_file} is a pedigree: founders with +both parents unknown are skipped instead of searched. Default is NULL.} } \value{ A named list (returned invisibly) with elements: @@ -123,6 +106,23 @@ Assigns the most likely parent(s) to each progeny from SNP genotype data using Mendelian error rates or homozygous mismatch rates. Parents or progeny absent from the genotype file are removed with a warning. } +\details{ +\strong{Method "fill_pedigree" ("fill in the blanks").} The +\code{progeny_file} must have \code{male_parent} and \code{female_parent} columns. A parent coded \code{0}, \code{NA} or +empty is treated as unknown and each row is handled on its own: +a known male searches only females, a known female searches only males, +and a row with both unknown gets the full pair search (unless listed in +\code{founders_file}). Rows with both parents known are skipped. The search is +the "best_pair" Mendelian trio test with the known parent held +fixed, and candidates come from \code{parents_file} (sex M/F/A as before); the +known parent only needs genotypes. Two columns are added to the results: +\code{search_mode} ("pair", "known_male" or "known_female") and +\code{known_parent_error_pct} (homozygous mismatch between the known parent and +the progeny; a high value suggests the known parent is wrong). A known +parent absent from the genotype file gives +\code{status = "known_parent_no_genotype"}. Without \code{method = "fill_pedigree"}, +any \code{male_parent} / \code{female_parent} columns in \code{progeny_file} are ignored. +} \examples{ \donttest{ geno_df <- data.frame( diff --git a/man/validate_pedigree.Rd b/man/validate_pedigree.Rd index 35d5bcf..7fd114a 100644 --- a/man/validate_pedigree.Rd +++ b/man/validate_pedigree.Rd @@ -13,7 +13,8 @@ validate_pedigree( single_parent_error_threshold = 2, verbose = TRUE, plot_results = TRUE, - ploidy = 2 + ploidy = 2, + marker_trio_table = FALSE ) } \arguments{ @@ -54,6 +55,10 @@ test (assumes autopolyploid inheritance; conservative for allopolyploids). Odd ploidy (e.g. triploid), where balanced gametes are undefined, falls back to a model-free opposite-homozygote exclusion evaluated on homozygous-informative markers only (reduced power). Default is 2.} + +\item{marker_trio_table}{Logical. If TRUE, the returned list includes +\code{marker_trio_table}, a marker-by-trio table of per-marker statuses (see +Details). Default is FALSE.} } \value{ An invisible named list with the following elements: @@ -71,6 +76,10 @@ with both parents known and genotyped: \code{marker}, \code{trios_tested}, \code{trios_mismatch} and \code{mismatch_pct}, sorted worst first. Markers with a much higher rate than the rest (e.g. paralogous or miscalled markers) are candidates for removal before rerunning.} +\item{marker_trio_table}{If \code{marker_trio_table = TRUE}, a data.table with +one row per marker (\code{marker}) and one column per trio with both parents +known and genotyped, named \code{id_male_parent_female_parent}; otherwise NULL. +See Details for the status values.} \item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} } } @@ -83,6 +92,16 @@ supplied. Trios absent from the genotype file are retained as no_genotype_data. } \details{ +The cells of \code{marker_trio_table} take one of these values: \code{match}; +\code{mismatch_male} or \code{mismatch_female} (the offspring is impossible with that +parent alone, whatever the other parent contributes); \code{mismatch_both} (the +offspring is impossible with either parent alone); \code{mismatch_combination} +(each parent is compatible alone but the pair is not, for example a 0 x 0 +cross with a heterozygous offspring); \code{untestable} (odd ploidy, marker not +informative); \verb{missing_*} (no call, named after every individual without a +call: \code{missing_progeny}, \code{missing_male}, \code{missing_female}, \code{missing_male_female}, +\code{missing_progeny_male}, \code{missing_progeny_female} or \code{missing_progeny_male_female}). + The \code{recommended_correction} column of \code{full_results} takes these values: \describe{ diff --git a/tests/testthat/test-assign_composition_line.R b/tests/testthat/test-assign_composition_line.R new file mode 100644 index 0000000..03dab6e --- /dev/null +++ b/tests/testthat/test-assign_composition_line.R @@ -0,0 +1,42 @@ +# tests/testthat/test-assign_composition_line.R + +comp <- rbind( + S1 = c(A = 0.998, B = 0.000, C = 0.002, R2 = 0.99), + S2 = c(A = 0.529, B = 0.471, C = 0.000, R2 = 0.97), + S3 = c(A = 0.340, B = 0.330, C = 0.330, R2 = 0.95) +) + +test_that("default threshold is 50 and weak assignments are Undetermined", { + out <- assign_composition_line(comp) + expect_equal(names(out), c("ID", "A", "B", "C", "Predicted line")) + expect_equal(out$ID, c("S1", "S2", "S3")) + expect_equal(out$`Predicted line`, c("A", "A", "Undetermined")) + expect_equal(out$A, c(99.8, 52.9, 34.0)) +}) + +test_that("threshold 0 always assigns the highest population", { + out <- assign_composition_line(comp, assign_threshold = 0) + expect_equal(out$`Predicted line`, c("A", "A", "A")) +}) + +test_that("a proportion equal to the threshold is assigned", { + m <- rbind(S1 = c(A = 0.5, B = 0.5)) + expect_equal(assign_composition_line(m, 50)$`Predicted line`, "A") +}) + +test_that("R2 is dropped and invalid thresholds are rejected", { + expect_false("R2" %in% names(assign_composition_line(comp))) + expect_error(assign_composition_line(comp, assign_threshold = 101)) + expect_error(assign_composition_line(comp, assign_threshold = "a")) +}) + +test_that("works on solve_composition_poly() output", { + X <- matrix(c(0.625, 0.5, 0.5, 0.5, 0.5, 0.5, 0.75, 0.5, 0.625, 0.625), + nrow = 5, byrow = TRUE, + dimnames = list(paste0("SNP", 1:5), c("VarA", "VarB"))) + Y <- matrix(c(2, 1, 2, 3, 4, 3, 4, 2, 3, 0), nrow = 2, byrow = TRUE, + dimnames = list(paste0("Test", 1:2), paste0("SNP", 1:5))) + out <- assign_composition_line(solve_composition_poly(Y, X, ploidy = 4)) + expect_equal(nrow(out), 2) + expect_true(all(out$`Predicted line` %in% c("VarA", "VarB", "Undetermined"))) +}) diff --git a/tests/testthat/test-validate_pedigree.R b/tests/testthat/test-validate_pedigree.R index 38c6288..ba31a9c 100644 --- a/tests/testthat/test-validate_pedigree.R +++ b/tests/testthat/test-validate_pedigree.R @@ -151,7 +151,7 @@ test_that("returns an invisible named list with all required elements", { expect_type(out, "list") expect_named(out, c("pass", "fail", "low_markers", "no_genotype_data", "founders", "missing_parents", "full_results", - "corrected_pedigree", "plot"), + "corrected_pedigree", "marker_summary", "marker_trio_table", "plot"), ignore.order = TRUE) }) @@ -605,7 +605,7 @@ test_that("verbose = TRUE returns valid named list without error", { expect_type(out, "list") expect_named(out, c("pass", "fail", "low_markers", "no_genotype_data", "founders", "missing_parents", "full_results", - "corrected_pedigree", "plot"), + "corrected_pedigree", "marker_summary", "marker_trio_table", "plot"), ignore.order = TRUE) }) @@ -735,3 +735,31 @@ test_that("invalid input type raises an error for validate_pedigree", { regexp = "Error reading input files" ) }) + +# ============================================================================== +# low_markers_* decisions remove only the identified parent +# ============================================================================== +test_that("corrected_pedigree: low_markers_remove_female_parent keeps the male parent", { + # IND_E (all 0) x IND_A (all 0) is fine for the male parent, but the female + # parent IND_B (all 2) is incompatible; only 4 markers -> low_markers + genos <- make_genos()[, 1:5] + ped <- data.table(id = "IND_E", male_parent = "IND_A", female_parent = "IND_B") + out <- validate_pedigree(ped, genos, min_markers = 10, + verbose = FALSE, plot_results = FALSE) + r <- out$full_results[id == "IND_E"] + expect_equal(r$status, "low_markers") + expect_equal(r$recommended_correction, "low_markers_remove_female_parent") + expect_equal(out$corrected_pedigree[id == "IND_E"]$female_parent, "0") + expect_equal(out$corrected_pedigree[id == "IND_E"]$male_parent, "IND_A") +}) + +test_that("corrected_pedigree: low_markers_remove_male_parent keeps the female parent", { + genos <- make_genos()[, 1:5] + ped <- data.table(id = "IND_E", male_parent = "IND_B", female_parent = "IND_A") + out <- validate_pedigree(ped, genos, min_markers = 10, + verbose = FALSE, plot_results = FALSE) + r <- out$full_results[id == "IND_E"] + expect_equal(r$recommended_correction, "low_markers_remove_male_parent") + expect_equal(out$corrected_pedigree[id == "IND_E"]$male_parent, "0") + expect_equal(out$corrected_pedigree[id == "IND_E"]$female_parent, "IND_A") +}) diff --git a/tests/testthat/test-validate_pedigree_marker_trio.R b/tests/testthat/test-validate_pedigree_marker_trio.R new file mode 100644 index 0000000..3fb8df1 --- /dev/null +++ b/tests/testthat/test-validate_pedigree_marker_trio.R @@ -0,0 +1,55 @@ +# tests/testthat/test-validate_pedigree_marker_trio.R + +library(testthat) +library(data.table) + +mt_data <- function() { + # one trio (C = P1 x P2), 11 markers, one status per marker + geno <- data.frame( + id = c("P1", "P2", "C"), + M1 = c(0L, 2L, 1L), # match + M2 = c(0L, 2L, 0L), # offspring impossible with the female alone + M3 = c(0L, 0L, 2L), # impossible with either parent + M4 = c(2L, 0L, 0L), # impossible with the male alone + M5 = c(NA, 1L, 1L), # male missing + M6 = c(0L, 2L, NA), # progeny missing + M7 = c(0L, NA, 1L), # female missing + M8 = c(0L, 0L, 1L), # each parent fine alone, pair impossible + M9 = c(NA, NA, 1L), # both parents missing + M10 = c(NA, 1L, NA), # progeny and male missing + M11 = c(NA, NA, NA) # all missing + ) + ped <- data.frame(id = "C", male_parent = "P1", female_parent = "P2", + stringsAsFactors = FALSE) + list(geno = geno, ped = ped) +} + +test_that("marker_trio_table is NULL unless requested", { + d <- mt_data() + out <- validate_pedigree(d$ped, d$geno, verbose = FALSE, plot_results = FALSE) + expect_null(out$marker_trio_table) +}) + +test_that("marker_trio_table gives the expected status per marker", { + d <- mt_data() + out <- validate_pedigree(d$ped, d$geno, verbose = FALSE, plot_results = FALSE, + marker_trio_table = TRUE) + tab <- as.data.frame(out$marker_trio_table) + expect_equal(names(tab), c("marker", "C_P1_P2")) + expect_equal(tab$marker, paste0("M", 1:11)) + expect_equal(tab$C_P1_P2, c( + "match", "mismatch_female", "mismatch_both", "mismatch_male", + "missing_male", "missing_progeny", "missing_female", + "mismatch_combination", "missing_male_female", + "missing_progeny_male", "missing_progeny_male_female")) +}) + +test_that("one column per trio with both parents known and genotyped", { + d <- mt_data() + d$ped <- rbind(d$ped, + data.frame(id = "P1", male_parent = "0", female_parent = "0"), + data.frame(id = "X", male_parent = "P1", female_parent = "P2")) + out <- validate_pedigree(d$ped, d$geno, verbose = FALSE, plot_results = FALSE, + marker_trio_table = TRUE) + expect_equal(setdiff(names(out$marker_trio_table), "marker"), "C_P1_P2") +}) diff --git a/vignettes/BIGpopA.Rmd b/vignettes/BIGpopA.Rmd index cbe3785..da847ff 100644 --- a/vignettes/BIGpopA.Rmd +++ b/vignettes/BIGpopA.Rmd @@ -280,21 +280,14 @@ prediction It is important to note that while the model produces an R2 value, it is expected to be low since an allele dosage (discrete number) is being regressed on an allele frequency (continuous value). -Assign each sample to the cultivar with the highest composition (ignoring the -R2 fit column) and format the panel columns as percentages: - +`assign_composition_line()` turns the result into the table used by Familia: +the estimated percentage of each panel column 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`. Use +`assign_threshold = 0` to always assign the highest line. ```{r} -line_cols = setdiff(names(prediction), "R2") - -pred_results = prediction %>% - mutate( - ID = rownames(prediction), - predicted_line = line_cols[max.col(as.matrix(prediction[line_cols]), ties.method = "first")], - across(all_of(line_cols), ~ scales::percent(., accuracy = 0.1)) - ) %>% - select(ID, predicted_line, all_of(line_cols)) - +pred_results = assign_composition_line(solve_composition_poly(validation, freq, ploidy = 4)) pred_results ``` @@ -314,6 +307,7 @@ assigned = find_parentage(genotypes, parents, progeny, method = "best_pair")$fu # breed/line composition 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) ``` # Session information From 4bd76b73193fe110e851853cb44a1070c1ac72f9 Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Mon, 12 Oct 2026 10:03:26 -0400 Subject: [PATCH 5/8] added option for all inputs to be piped internally --- NEWS.md | 1 + R/find_parentage.R | 9 +++------ R/utils.R | 18 ++++++++++++++++++ R/validate_pedigree.R | 15 +++------------ _to_delete/test-founders_input.R | 9 +++++++++ man/find_parentage.Rd | 4 ++-- man/validate_pedigree.Rd | 4 ++-- tests/testthat/test-founders_input.R | 9 +++++++++ 8 files changed, 47 insertions(+), 22 deletions(-) create mode 100644 _to_delete/test-founders_input.R create mode 100644 tests/testthat/test-founders_input.R diff --git a/NEWS.md b/NEWS.md index 66f583c..56ae501 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,5 @@ # BIGpopA 2.2.0 +- `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. - `find_parentage()`: new `method = "fill_pedigree"` fills in missing parents of a pedigree. `progeny_file` must have `id`, `male_parent` and `female_parent` columns (e.g. `validate_pedigree()$corrected_pedigree`); parents coded `0`/`NA`/empty are unknown. A known male searches only females, a known female only males, and rows with both unknown get the full pair search (skip founders with the new `founders_file`). Rows with both parents known are skipped. Results gain `search_mode` and `known_parent_error_pct`, and a known parent without genotypes gets `status = "known_parent_no_genotype"`. diff --git a/R/find_parentage.R b/R/find_parentage.R index e679ff9..216a96b 100644 --- a/R/find_parentage.R +++ b/R/find_parentage.R @@ -42,8 +42,8 @@ #' Odd ploidy (e.g. triploid), where balanced gametes are undefined, falls #' back to a model-free opposite-homozygote exclusion evaluated on #' homozygous-informative markers only (reduced power). Default is 2. -#' @param founders_file Character, optional. Path to a one-column file of -#' founder IDs. Only used when `progeny_file` is a pedigree: founders with +#' @param founders_file Optional. Founder IDs as a path to a one-column file, a +#' character vector of ids, or a data.frame (first column). Only used when `progeny_file` is a pedigree: founders with #' both parents unknown are skipped instead of searched. Default is NULL. #' @param top_n Integer >= 1. Number of best candidates to report per progeny. #' With the default of 1, only the best assignment is reported (plus tied @@ -199,10 +199,7 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, method <- "best_pair" progeny_candidates[, male_parent := base::as.character(male_parent)] progeny_candidates[, female_parent := base::as.character(female_parent)] - founder_ids <- if (!is.null(founders_file)) - base::unique(data.table::fread(founders_file, header = FALSE, - colClasses = "character")[[1]]) else - base::character(0) + founder_ids <- .read_founders(founders_file) um <- is_unknown(progeny_candidates$male_parent) uf <- is_unknown(progeny_candidates$female_parent) is_founder <- um & uf & progeny_candidates$id %in% founder_ids diff --git a/R/utils.R b/R/utils.R index 6027216..c5efe63 100644 --- a/R/utils.R +++ b/R/utils.R @@ -349,3 +349,21 @@ QPsolve <- function(Y, X) { TRUE ~ "mismatch_combination" ) } + + +# Read founder ids from a file path, a character vector of ids, or a +# data.frame / data.table (first column). Returns a character vector. +.read_founders <- function(x) { + if (is.null(x)) return(character(0)) + if (is.data.frame(x)) return(unique(as.character(x[[1]]))) + if (is.character(x) && length(x) == 1 && file.exists(x)) { + out <- tryCatch( + data.table::fread(x, header = FALSE, colClasses = "character"), + error = function(e) + stop("Could not read founders list. Ensure it is a plain text or CSV/TSV file.") + ) + return(unique(out[[1]])) + } + if (is.atomic(x)) return(unique(as.character(x))) + stop("founders_file must be a file path, a vector of ids, or a data.frame.") +} diff --git a/R/validate_pedigree.R b/R/validate_pedigree.R index ea723fb..00ce2f5 100644 --- a/R/validate_pedigree.R +++ b/R/validate_pedigree.R @@ -15,8 +15,8 @@ #' (diploid only), converted with [ped_to_dosage()]; or a data.frame / #' data.table with an id column followed by marker columns coded as #' allele-B dosage (0, 1, ..., ploidy). -#' @param founders_file Character, optional. Path to a one-column file listing -#' founder IDs. Founders with both parents coded as 0 are left unchanged. +#' @param founders_file Optional. Founder IDs as a path to a one-column file, a +#' character vector of ids, or a data.frame (first column). Founders with both parents coded as 0 are left unchanged. #' Defaults to NULL. #' @param trio_error_threshold Numeric. Maximum Mendelian error percentage to #' classify a trio as pass (default: 5.0). Must be between 0 and 100. @@ -190,16 +190,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, original_pedigree <- data.table::copy(pedigree) #### Read founders list #### - if (!is.null(founders_file)) { - founders_raw <- tryCatch({ - data.table::fread(founders_file, header = FALSE, colClasses = "character") - }, error = function(e) { - stop("Could not read founders list. Ensure it is a plain text or CSV/TSV file.") - }) - founder_ids <- unique(founders_raw[[1]]) - } else { - founder_ids <- character(0) - } + founder_ids <- .read_founders(founders_file) #### Build genotype matrices #### genos_mat <- base::as.matrix(genos, rownames = "id") diff --git a/_to_delete/test-founders_input.R b/_to_delete/test-founders_input.R new file mode 100644 index 0000000..b46da85 --- /dev/null +++ b/_to_delete/test-founders_input.R @@ -0,0 +1,9 @@ +test_that(".read_founders accepts NULL, vector, data.frame and file", { + expect_equal(BIGpopA:::.read_founders(NULL), character(0)) + expect_equal(BIGpopA:::.read_founders(c("A", "B", "A")), c("A", "B")) + expect_equal(BIGpopA:::.read_founders(data.frame(id = c("A", "B"))), c("A", "B")) + f <- tempfile(fileext = ".txt") + writeLines(c("id", "A", "B"), f) + expect_equal(BIGpopA:::.read_founders(f), c("A", "B")) + expect_error(BIGpopA:::.read_founders(list(1)), "founders_file") +}) diff --git a/man/find_parentage.Rd b/man/find_parentage.Rd index b8b9b27..08e1eae 100644 --- a/man/find_parentage.Rd +++ b/man/find_parentage.Rd @@ -84,8 +84,8 @@ for the other methods). Ties are then ranked like any other candidate, so \code{show_ties} has no effect. Progeny with fewer than \code{top_n} candidates have \code{NA} in the remaining ranks. Status refers to rank 1 only.} -\item{founders_file}{Character, optional. Path to a one-column file of -founder IDs. Only used when \code{progeny_file} is a pedigree: founders with +\item{founders_file}{Optional. Founder IDs as a path to a one-column file, a +character vector of ids, or a data.frame (first column). Only used when \code{progeny_file} is a pedigree: founders with both parents unknown are skipped instead of searched. Default is NULL.} } \value{ diff --git a/man/validate_pedigree.Rd b/man/validate_pedigree.Rd index 7fd114a..1d0b2bc 100644 --- a/man/validate_pedigree.Rd +++ b/man/validate_pedigree.Rd @@ -28,8 +28,8 @@ with \code{\link[=vcf_to_dosage]{vcf_to_dosage()}} using \code{ploidy}; a path t data.table with an id column followed by marker columns coded as allele-B dosage (0, 1, ..., ploidy).} -\item{founders_file}{Character, optional. Path to a one-column file listing -founder IDs. Founders with both parents coded as 0 are left unchanged. +\item{founders_file}{Optional. Founder IDs as a path to a one-column file, a +character vector of ids, or a data.frame (first column). Founders with both parents coded as 0 are left unchanged. Defaults to NULL.} \item{trio_error_threshold}{Numeric. Maximum Mendelian error percentage to diff --git a/tests/testthat/test-founders_input.R b/tests/testthat/test-founders_input.R new file mode 100644 index 0000000..4003bea --- /dev/null +++ b/tests/testthat/test-founders_input.R @@ -0,0 +1,9 @@ +test_that(".read_founders accepts NULL, vector, data.frame and file", { + expect_equal(BIGpopA:::.read_founders(NULL), character(0)) + expect_equal(BIGpopA:::.read_founders(c("A", "B", "A")), c("A", "B")) + expect_equal(BIGpopA:::.read_founders(data.frame(id = c("A", "B"))), c("A", "B")) + f <- tempfile(fileext = ".txt") + writeLines(c("A", "B"), f) + expect_equal(BIGpopA:::.read_founders(f), c("A", "B")) + expect_error(BIGpopA:::.read_founders(list(1)), "founders_file") +}) From fc41f11b14d5444a94f4d167f7f1f0a75f0ba59c Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Mon, 12 Oct 2026 10:14:07 -0400 Subject: [PATCH 6/8] added a set of hard tests made by Claude --- tests/testthat/test-hard.R | 371 +++++++++++++++++++++++++++++++++++++ 1 file changed, 371 insertions(+) create mode 100644 tests/testthat/test-hard.R diff --git a/tests/testthat/test-hard.R b/tests/testthat/test-hard.R new file mode 100644 index 0000000..7fa1217 --- /dev/null +++ b/tests/testthat/test-hard.R @@ -0,0 +1,371 @@ +# tests/testthat/test-hard.R +# Randomised, oracle-based and invariance tests. The oracles below are written +# independently of the package code (explicit gamete enumeration, brute force +# over all candidate pairs), so a pass means the package agrees with the +# definition of the test, not with itself. + +library(testthat) +library(data.table) + +# ============================================================================== +# Helpers +# ============================================================================== + +# gamete allele-B counts a parent with dosage g can transmit at ploidy P, +# found by explicit enumeration of k alleles drawn from g B and (P - g) A +oracle_gametes <- function(g, P) { + h <- P / 2 + Filter(function(k) k <= g && (h - k) <= (P - g), 0:h) +} + +# TRUE if offspring dosage o is possible from parent dosages gm and gf +oracle_valid <- function(o, gm, gf, P) { + km <- oracle_gametes(gm, P) + kf <- oracle_gametes(gf, P) + any(outer(km, kf, "+") == o) +} + +# returns c(tested, mismatches) for one trio over all markers (NA = skip) +oracle_trio <- function(o, m, f, P) { + ok <- !is.na(o) & !is.na(m) & !is.na(f) + bad <- vapply(which(ok), function(j) !oracle_valid(o[j], m[j], f[j], P), logical(1)) + c(tested = sum(ok), mismatches = sum(bad)) +} + +# simulate a child dosage from two parents (random meiosis) +sim_child <- function(m, f, P) { + gam <- function(g) sum(sample(rep(c(1L, 0L), c(g, P - g)), P / 2)) + vapply(seq_along(m), function(j) gam(m[j]) + gam(f[j]), numeric(1)) +} + +# random panel of parents, valid and random children, with optional NA +make_panel <- function(P, n_markers = 60, n_trios = 24, na_rate = 0.1, seed = 1) { + set.seed(seed) + mk <- paste0("S", seq_len(n_markers)) + par_ids <- paste0("P", seq_len(2 * n_trios)) + par <- matrix(sample(0:P, length(par_ids) * n_markers, replace = TRUE), + nrow = length(par_ids), dimnames = list(par_ids, mk)) + kids <- matrix(NA_real_, nrow = n_trios, ncol = n_markers, + dimnames = list(paste0("C", seq_len(n_trios)), mk)) + ped <- data.frame(id = rownames(kids), + male_parent = par_ids[seq(1, 2 * n_trios, by = 2)], + female_parent = par_ids[seq(2, 2 * n_trios, by = 2)], + stringsAsFactors = FALSE) + for (i in seq_len(n_trios)) { + m <- par[ped$male_parent[i], ] + f <- par[ped$female_parent[i], ] + kids[i, ] <- if (i %% 2 == 0) sim_child(m, f, P) # valid trio + else sample(0:P, n_markers, replace = TRUE) # random child + } + all_g <- rbind(par, kids) + if (na_rate > 0) + all_g[sample(length(all_g), round(na_rate * length(all_g)))] <- NA + geno <- data.frame(id = rownames(all_g), all_g, check.names = FALSE, + stringsAsFactors = FALSE, row.names = NULL) + list(geno = geno, ped = ped, mat = all_g) +} + +quiet_validate <- function(ped, geno, ...) { + validate_pedigree(ped, geno, verbose = FALSE, plot_results = FALSE, ...) +} + +# ============================================================================== +# 1. validate_pedigree agrees with brute-force Mendelian oracle +# ============================================================================== + +for (P in c(2, 4, 6, 8)) { + test_that(sprintf("validate_pedigree matches the gamete oracle (ploidy %d, with NA)", P), { + d <- make_panel(P, seed = 100 + P) + out <- quiet_validate(d$ped, d$geno, ploidy = P, min_markers = 3)$full_results + out <- as.data.frame(out) + + for (i in seq_len(nrow(d$ped))) { + o <- oracle_trio(d$mat[d$ped$id[i], ], + d$mat[d$ped$male_parent[i], ], + d$mat[d$ped$female_parent[i], ], P) + r <- out[out$id == d$ped$id[i], ] + expect_equal(r$trio_markers_tested, unname(o["tested"])) + expect_equal(r$trio_mendelian_error_pct, + round(100 * unname(o["mismatches"]) / unname(o["tested"]), 2), + tolerance = 0.011) + expect_equal(r$status, + if (100 * o["mismatches"] / o["tested"] <= 5) "pass" else "fail") + } + }) +} + +test_that("valid simulated trios never fail (no false positives, no NA)", { + for (P in c(2, 4, 6, 8)) { + d <- make_panel(P, na_rate = 0, seed = 7 + P) + out <- as.data.frame(quiet_validate(d$ped, d$geno, ploidy = P)$full_results) + valid_ids <- d$ped$id[seq(2, nrow(d$ped), by = 2)] + expect_true(all(out$trio_mendelian_error_pct[out$id %in% valid_ids] == 0), + info = paste("ploidy", P)) + } +}) + +# ============================================================================== +# 2. Odd ploidy (hand-computed) +# ============================================================================== + +test_that("triploid trio uses homozygous opposite-allele test only", { + geno <- data.frame( + id = c("M", "F", "C"), + M1 = c(3L, 0L, 0L), # male 3 vs offspring 0 -> testable, mismatch + M2 = c(0L, 0L, 0L), # testable, match + M3 = c(3L, 1L, 3L), # male homozygous same as offspring -> match + M4 = c(1L, 2L, 1L) # offspring heterozygous -> untestable + ) + ped <- data.frame(id = "C", male_parent = "M", female_parent = "F") + r <- as.data.frame(quiet_validate(ped, geno, ploidy = 3, min_markers = 1)$full_results) + expect_equal(r$trio_markers_tested, 3) + expect_equal(r$trio_mendelian_error_pct, round(100 / 3, 2), tolerance = 0.011) +}) + +# ============================================================================== +# 3. Invariance: marker order, trio order and input type do not matter +# ============================================================================== + +test_that("results do not depend on marker order or trio order", { + d <- make_panel(4, seed = 33) + base <- as.data.frame(quiet_validate(d$ped, d$geno, ploidy = 4)$full_results) + set.seed(5) + geno2 <- d$geno[sample(nrow(d$geno)), c("id", sample(setdiff(names(d$geno), "id")))] + ped2 <- d$ped[sample(nrow(d$ped)), ] + perm <- as.data.frame(quiet_validate(ped2, geno2, ploidy = 4)$full_results) + cols <- c("id", "trio_mendelian_error_pct", "trio_markers_tested", "status") + expect_equal(base[order(base$id), cols], perm[order(perm$id), cols], + ignore_attr = TRUE) +}) + +test_that("file input and in-memory input give identical results", { + d <- make_panel(4, seed = 34) + fg <- tempfile(fileext = ".txt"); fp <- tempfile(fileext = ".txt") + fwrite(d$geno, fg, sep = "\t"); fwrite(d$ped, fp, sep = "\t") + a <- as.data.frame(quiet_validate(d$ped, d$geno, ploidy = 4)$full_results) + b <- as.data.frame(quiet_validate(fp, fg, ploidy = 4)$full_results) + expect_equal(a, b, ignore_attr = TRUE) +}) + +test_that("founders as file, vector and data.frame give the same result", { + d <- make_panel(2, seed = 35) + ped <- rbind(d$ped, data.frame(id = c("P1", "P2"), male_parent = "0", female_parent = "0")) + ff <- tempfile(fileext = ".txt"); writeLines(c("P1", "P2"), ff) + a <- quiet_validate(ped, d$geno, founders_file = ff)$full_results + b <- quiet_validate(ped, d$geno, founders_file = c("P1", "P2"))$full_results + cc <- quiet_validate(ped, d$geno, founders_file = data.frame(id = c("P1", "P2")))$full_results + expect_equal(a, b) + expect_equal(a, cc) + expect_equal(sum(a$status == "founders"), 2) +}) + +# ============================================================================== +# 4. marker_trio_table consistent with the oracle +# ============================================================================== + +test_that("marker_trio_table mismatch and missing counts agree with the oracle", { + P <- 4 + d <- make_panel(P, seed = 41) + out <- quiet_validate(d$ped, d$geno, ploidy = P, marker_trio_table = TRUE) + tab <- as.data.frame(out$marker_trio_table) + expect_equal(tab$marker, setdiff(names(d$geno), "id")) + + for (i in seq_len(nrow(d$ped))) { + col <- paste(d$ped$id[i], d$ped$male_parent[i], d$ped$female_parent[i], sep = "_") + lab <- tab[[col]] + o <- d$mat[d$ped$id[i], ] + m <- d$mat[d$ped$male_parent[i], ] + f <- d$mat[d$ped$female_parent[i], ] + any_na <- unname(is.na(o) | is.na(m) | is.na(f)) + expect_equal(grepl("^missing", lab), any_na) + expect_equal(sum(grepl("^mismatch", lab)), + unname(oracle_trio(o, m, f, P)["mismatches"])) + expect_equal(sum(lab == "match"), + unname(oracle_trio(o, m, f, P)["tested"] - + oracle_trio(o, m, f, P)["mismatches"])) + } +}) + +# ============================================================================== +# 5. marker_summary consistent with the trio results +# ============================================================================== + +test_that("marker_summary totals equal the per-trio totals", { + d <- make_panel(4, seed = 51) + out <- quiet_validate(d$ped, d$geno, ploidy = 4) + ms <- as.data.frame(out$marker_summary) + fr <- as.data.frame(out$full_results) + tested_trios <- fr[fr$trio_markers_tested > 0, ] + expect_equal(sum(ms$trios_tested), sum(tested_trios$trio_markers_tested)) + expect_equal(sum(ms$trios_mismatch), + sum(round(tested_trios$trio_markers_tested * + tested_trios$trio_mendelian_error_pct / 100))) + expect_true(all(ms$trios_mismatch <= ms$trios_tested)) + expect_false(is.unsorted(rev(ms$mismatch_pct[!is.na(ms$mismatch_pct)]))) +}) + +# ============================================================================== +# 6. find_parentage best_pair against brute force +# ============================================================================== + +test_that("best_pair error equals the brute-force minimum over all pairs", { + skip_on_cran() + P <- 4 + set.seed(61) + n_m <- 50; n_male <- 5; n_fem <- 5 + mk <- paste0("S", seq_len(n_m)) + males <- matrix(sample(0:P, n_male * n_m, TRUE), n_male, + dimnames = list(paste0("MA", 1:n_male), mk)) + fems <- matrix(sample(0:P, n_fem * n_m, TRUE), n_fem, + dimnames = list(paste0("FE", 1:n_fem), mk)) + true_pairs <- data.frame(id = paste0("K", 1:6), + m = rep(rownames(males), length.out = 6)[c(1, 2, 3, 4, 5, 1)], + f = rep(rownames(fems), length.out = 6)[c(2, 3, 4, 5, 1, 4)], + stringsAsFactors = FALSE) + kids <- t(vapply(seq_len(nrow(true_pairs)), function(i) + sim_child(males[true_pairs$m[i], ], fems[true_pairs$f[i], ], P), numeric(n_m))) + dimnames(kids) <- list(true_pairs$id, mk) + all_g <- rbind(males, fems, kids) + all_g[sample(length(all_g), 15)] <- NA # a few missing calls + geno <- data.frame(id = rownames(all_g), all_g, check.names = FALSE, row.names = NULL) + parents <- data.frame(id = c(rownames(males), rownames(fems)), + sex = rep(c("M", "F"), c(n_male, n_fem))) + progeny <- data.frame(id = true_pairs$id) + + res <- as.data.frame(find_parentage(geno, parents, progeny, method = "best_pair", + ploidy = P, min_markers = 3, verbose = FALSE, + plot_results = FALSE)$full_results) + + for (k in true_pairs$id) { + errs <- outer(rownames(males), rownames(fems), Vectorize(function(a, b) { + o <- oracle_trio(all_g[k, ], all_g[a, ], all_g[b, ], P) + 100 * o["mismatches"] / o["tested"] + })) + expect_equal(as.numeric(res$mendelian_error_pct[res$id == k]), round(min(errs), 2), + tolerance = 0.011, info = k) + expect_equal(res$tied_candidates[res$id == k], sum(abs(errs - min(errs)) < 1e-9), + info = k) + } + # the simulated pair always has the lowest error unless another pair ties + unique_best <- res$tied_candidates == 1 + expect_equal(res$male_parent[unique_best], true_pairs$m[match(res$id[unique_best], true_pairs$id)]) + expect_equal(res$female_parent[unique_best], true_pairs$f[match(res$id[unique_best], true_pairs$id)]) +}) + +test_that("top_n = 3 returns non-decreasing errors per progeny", { + P <- 4 + d <- make_panel(P, na_rate = 0, seed = 71) + pid <- grep("^P", d$geno$id, value = TRUE) + parents <- data.frame(id = pid, sex = rep(c("M", "F"), length.out = length(pid))) + progeny <- data.frame(id = d$ped$id[1:4]) + r <- as.data.frame(find_parentage(d$geno, parents, progeny, method = "best_pair", + ploidy = P, top_n = 3, verbose = FALSE, + plot_results = FALSE)$full_results) + e <- r[, c("mendelian_error_pct", "mendelian_error_pct_2", "mendelian_error_pct_3")] + e[] <- lapply(e, as.numeric) + expect_true(all(e[, 1] <= e[, 2] & e[, 2] <= e[, 3])) +}) + +# ============================================================================== +# 7. check_ped invariants on a larger random pedigree +# ============================================================================== + +make_big_ped <- function(n = 400, seed = 81) { + set.seed(seed) + ids <- sprintf("I%04d", seq_len(n)) + sires <- ids[seq(1, n, by = 2)] # fixed sex roles + dams <- ids[seq(2, n, by = 2)] + ped <- data.frame(id = ids, male_parent = "0", female_parent = "0", + stringsAsFactors = FALSE) + for (i in 21:n) { # first 20 are founders; parents come earlier + ped$male_parent[i] <- sample(sires[sires %in% ids[seq_len(i - 1)]], 1) + ped$female_parent[i] <- sample(dams[dams %in% ids[seq_len(i - 1)]], 1) + } + ped +} + +test_that("check_ped on a clean random pedigree changes nothing", { + ped <- make_big_ped() + out <- check_ped(ped, verbose = FALSE)$corrected_pedigree + key <- function(x) sort(paste(x$id, x$male_parent, x$female_parent)) + expect_equal(key(out), key(ped)) +}) + +test_that("check_ped removes injected duplicates and is idempotent", { + ped <- make_big_ped() + dup <- ped[sample(nrow(ped), 40), ] + messy <- rbind(ped, dup)[sample(nrow(ped) + 40), ] + once <- check_ped(messy, verbose = FALSE)$corrected_pedigree + twice <- check_ped(once, verbose = FALSE)$corrected_pedigree + key <- function(x) sort(paste(x$id, x$male_parent, x$female_parent)) + expect_equal(nrow(once), nrow(ped)) + expect_false(any(duplicated(once))) + expect_equal(key(once), key(ped)) + expect_equal(key(twice), key(once)) +}) + +test_that("every non-zero parent has its own row after check_ped", { + ped <- make_big_ped() + ped$male_parent[ped$id == "I0100"] <- "I0003" # make sure I0003 is used as a parent + ped <- ped[ped$id != "I0003", ] # then drop its own row + out <- check_ped(ped, verbose = FALSE)$corrected_pedigree + parents <- setdiff(c(out$male_parent, out$female_parent), "0") + expect_true(all(parents %in% out$id)) + expect_true("I0003" %in% out$id) +}) + +test_that("a cycle inside a large pedigree is reported and nothing is dropped", { + ped <- make_big_ped() + ped$male_parent[ped$id == "I0399"] <- "I0001" # two sires that are each + ped$male_parent[ped$id == "I0001"] <- "I0399" # other's sire: a cycle + res <- check_ped(ped, verbose = FALSE) + expect_gt(nrow(res$dependencies), 0) + expect_true(all(ped$id %in% res$corrected_pedigree$id)) +}) + +# ============================================================================== +# 8. Composition: exact recovery of known mixtures +# ============================================================================== + +test_that("solve_composition_poly recovers noise-free mixtures", { + set.seed(91) + P <- 4 + nm <- 80 + # lines are fixed for one allele at every marker: line frequencies are 0 or 1 + a <- sample(0:1, nm, replace = TRUE) + X <- cbind(LineA = a, LineB = 1 - a) + rownames(X) <- paste0("S", seq_len(nm)) + w <- rbind(c(1, 0), c(0, 1), c(0.5, 0.5), c(0.25, 0.75), c(0.75, 0.25)) + Y <- t(apply(w, 1, function(wi) P * as.vector(X %*% wi))) # integer dosages + dimnames(Y) <- list(paste0("A", 1:5), rownames(X)) + est <- suppressWarnings(solve_composition_poly(Y, X, ploidy = P)) # R2 warns when sd = 0 + expect_equal(unname(est[, 1:2]), w, tolerance = 1e-3) + expect_equal(unname(rowSums(est[, 1:2])), rep(1, 5), tolerance = 1e-8) +}) + +test_that("allele_freq_poly matches colMeans / ploidy for random data", { + set.seed(92) + P <- 6 + g <- matrix(sample(0:P, 30 * 40, TRUE), 30, + dimnames = list(paste0("I", 1:30), paste0("S", 1:40))) + pops <- list(One = paste0("I", 1:10), Two = paste0("I", 11:30)) + f <- allele_freq_poly(g, pops, ploidy = P) + expect_equal(unname(f[, "One"]), unname(colMeans(g[1:10, ]) / P)) + expect_equal(unname(f[, "Two"]), unname(colMeans(g[11:30, ]) / P)) + expect_true(all(f >= 0 & f <= 1)) +}) + +# ============================================================================== +# 9. Scale +# ============================================================================== + +test_that("validate_pedigree handles 300 trios x 3000 markers (ploidy 4) in time", { + skip_on_cran() + d <- make_panel(4, n_markers = 3000, n_trios = 300, na_rate = 0.02, seed = 99) + t0 <- Sys.time() + out <- quiet_validate(d$ped, d$geno, ploidy = 4) + secs <- as.numeric(difftime(Sys.time(), t0, units = "secs")) + expect_equal(nrow(out$full_results), 300) + expect_equal(nrow(out$marker_summary), 3000) + expect_lt(secs, 120) +}) From 6588a5be52507248651d2f73dec308bbefee9155 Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Mon, 12 Oct 2026 10:20:27 -0400 Subject: [PATCH 7/8] deep code check by Claude --- R/allele_freq_poly.R | 2 +- R/find_parentage.R | 326 +++++++++++++++++++++--------------------- R/ped_to_dosage.R | 102 ++++++------- R/utils.R | 94 ++++++------ R/validate_pedigree.R | 149 +++++++++---------- R/vcf_to_dosage.R | 144 +++++++++---------- man/find_parentage.Rd | 2 +- man/ped_to_dosage.Rd | 2 +- 8 files changed, 411 insertions(+), 410 deletions(-) diff --git a/R/allele_freq_poly.R b/R/allele_freq_poly.R index a529396..b2bf114 100644 --- a/R/allele_freq_poly.R +++ b/R/allele_freq_poly.R @@ -51,7 +51,7 @@ allele_freq_poly <- function(geno, populations, ploidy = 2) { # Accept text files, VCFs, PLINK .ped, and in-memory tables (individuals x SNPs) geno <- .read_genotypes(geno, ploidy = ploidy, format = "matrix") - counted <- base::attr(geno, "counted_allele") # set only for .ped input + counted <- attr(geno, "counted_allele") # set only for .ped input # Initialize returned df X <- matrix(NA, nrow = ncol(geno), ncol = length(populations)) diff --git a/R/find_parentage.R b/R/find_parentage.R index 216a96b..fc09634 100644 --- a/R/find_parentage.R +++ b/R/find_parentage.R @@ -44,7 +44,7 @@ #' homozygous-informative markers only (reduced power). Default is 2. #' @param founders_file Optional. Founder IDs as a path to a one-column file, a #' character vector of ids, or a data.frame (first column). Only used when `progeny_file` is a pedigree: founders with -#' both parents unknown are skipped instead of searched. Default is NULL. +#' both parents unknown are skipped and not searched. Default is NULL. #' @param top_n Integer >= 1. Number of best candidates to report per progeny. #' With the default of 1, only the best assignment is reported (plus tied #' best pairs when `show_ties = TRUE` and `method = "best_pair"`). With @@ -156,9 +156,9 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, stop("error_threshold must be between 0 and 100.") .check_ploidy(ploidy) if (!is.numeric(top_n) || length(top_n) != 1 || is.na(top_n) || - top_n < 1 || top_n != base::round(top_n)) + top_n < 1 || top_n != round(top_n)) stop("top_n must be a single positive integer.") - top_n <- base::as.integer(top_n) + top_n <- as.integer(top_n) # Accept file path OR in-memory data.frame / data.table read_flex <- function(x, label, ...) { @@ -189,164 +189,163 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, c("id", "male_parent", "female_parent")) #### Pedigree ("fill in the blanks") mode #### - is_unknown <- function(x) base::is.na(x) | base::trimws(base::as.character(x)) %in% c("", "0", "NA") + is_unknown <- function(x) is.na(x) | trimws(as.character(x)) %in% c("", "0", "NA") method_reported <- method pedigree_mode <- method == "fill_pedigree" if (pedigree_mode) { - if (!base::all(c("male_parent", "female_parent") %in% base::names(progeny_candidates))) + if (!all(c("male_parent", "female_parent") %in% names(progeny_candidates))) stop("method = 'fill_pedigree' requires a progeny_file with id, male_parent ", "and female_parent columns.") method <- "best_pair" - progeny_candidates[, male_parent := base::as.character(male_parent)] - progeny_candidates[, female_parent := base::as.character(female_parent)] + progeny_candidates[, male_parent := as.character(male_parent)] + progeny_candidates[, female_parent := as.character(female_parent)] founder_ids <- .read_founders(founders_file) um <- is_unknown(progeny_candidates$male_parent) uf <- is_unknown(progeny_candidates$female_parent) is_founder <- um & uf & progeny_candidates$id %in% founder_ids if (verbose) - base::cat("Pedigree mode:", base::sum(!um & !uf), "complete trios skipped,", - base::sum(is_founder), "founders skipped,", - base::sum((um | uf) & !is_founder), "rows with a missing parent to search.\n") + cat(sprintf("Pedigree mode: %d complete trios skipped, %d founders skipped, %d rows with a missing parent to search.\n", + sum(!um & !uf), sum(is_founder), sum((um | uf) & !is_founder))) progeny_candidates <- progeny_candidates[(um | uf) & !is_founder] - if (base::nrow(progeny_candidates) == 0) + if (nrow(progeny_candidates) == 0) stop("No progeny with a missing parent to fill in.") } valid_ids <- genos$id - removed_parents <- base::setdiff(all_parents$id, valid_ids) - if (base::length(removed_parents) > 0) { + removed_parents <- setdiff(all_parents$id, valid_ids) + if (length(removed_parents) > 0) { warning("The following parent IDs were not in the genotype file and will not be analyzed: ", paste(removed_parents, collapse = ", "), call. = FALSE) all_parents <- all_parents[id %in% valid_ids] } - removed_progeny <- base::setdiff(progeny_candidates$id, valid_ids) - if (base::length(removed_progeny) > 0) { + removed_progeny <- setdiff(progeny_candidates$id, valid_ids) + if (length(removed_progeny) > 0) { warning("The following progeny IDs were not in the genotype file and will not be analyzed: ", paste(removed_progeny, collapse = ", "), call. = FALSE) progeny_candidates <- progeny_candidates[id %in% valid_ids] } - if (!"sex" %in% base::colnames(all_parents)) { + if (!"sex" %in% colnames(all_parents)) { warning("No 'sex' column in parents file. All parents treated as ambiguous ('A').") all_parents[, sex := "A"] } - all_parents[, sex := base::toupper(sex)] - all_parents <- base::unique(all_parents, by = c("id", "sex")) - male_parent_candidates <- base::unique(all_parents[sex %in% c("M", "A", "NA"), .SD], + all_parents[, sex := toupper(sex)] + all_parents <- unique(all_parents, by = c("id", "sex")) + male_parent_candidates <- unique(all_parents[sex %in% c("M", "A", "NA"), .SD], by = "id") - female_parent_candidates <- base::unique(all_parents[sex %in% c("F", "A", "NA")], + female_parent_candidates <- unique(all_parents[sex %in% c("F", "A", "NA")], by = "id") - if (base::nrow(male_parent_candidates) == 0 && + if (nrow(male_parent_candidates) == 0 && method %in% c("best_male_parent", "best_pair")) warning("No valid male parent candidates remain after filtering.", call. = FALSE) - if (base::nrow(female_parent_candidates) == 0 && + if (nrow(female_parent_candidates) == 0 && method %in% c("best_female_parent", "best_pair")) warning("No valid female parent candidates remain after filtering.", call. = FALSE) - if (base::nrow(progeny_candidates) == 0) + if (nrow(progeny_candidates) == 0) stop("No valid progeny candidates remain after filtering.") #### Pre-compute genotype matrices once #### - genos_mat <- base::as.matrix(genos, rownames = "id") + genos_mat <- as.matrix(genos, rownames = "id") genos_hom <- data.table::copy(genos) - marker_cols <- base::setdiff(base::names(genos_hom), "id") + marker_cols <- setdiff(names(genos_hom), "id") for (col in marker_cols) - genos_hom[base::get(col) > 0 & base::get(col) < ploidy, (col) := NA_integer_] - genos_hom_mat <- base::as.matrix(genos_hom, rownames = "id") + genos_hom[get(col) > 0 & get(col) < ploidy, (col) := NA_integer_] + genos_hom_mat <- as.matrix(genos_hom, rownames = "id") #### Status helper #### assign_status <- function(markers, error_pct) { - base::ifelse(markers < min_markers, "low_markers", - base::ifelse(error_pct > error_threshold, "high_error", "pass")) + ifelse(markers < min_markers, "low_markers", + ifelse(error_pct > error_threshold, "high_error", "pass")) } #### Logic for Homozygous Matching Methods #### if (method %in% c("best_male_parent", "best_female_parent", "best_match")) { - parent_ids <- base::switch(method, + parent_ids <- switch(method, "best_male_parent" = male_parent_candidates$id, "best_female_parent" = female_parent_candidates$id, - "best_match" = base::union(male_parent_candidates$id, + "best_match" = union(male_parent_candidates$id, female_parent_candidates$id)) - parent_genos <- genos_hom_mat[base::rownames(genos_hom_mat) %in% parent_ids, , drop = FALSE] - progeny_genos <- genos_hom_mat[base::rownames(genos_hom_mat) %in% progeny_candidates$id, , drop = FALSE] + parent_genos <- genos_hom_mat[rownames(genos_hom_mat) %in% parent_ids, , drop = FALSE] + progeny_genos <- genos_hom_mat[rownames(genos_hom_mat) %in% progeny_candidates$id, , drop = FALSE] - n_progeny <- base::nrow(progeny_genos) + n_progeny <- nrow(progeny_genos) results_dt <- data.table::data.table( - id = base::rownames(progeny_genos), + id = rownames(progeny_genos), best_match = NA_character_, mendelian_error_pct = NA_real_, markers_tested = NA_integer_, status = NA_character_, tied_candidates = NA_integer_ ) - extra_list <- base::vector("list", n_progeny) + extra_list <- vector("list", n_progeny) - for (i in base::seq_len(n_progeny)) { - prog_id <- base::rownames(progeny_genos)[i] + for (i in seq_len(n_progeny)) { + prog_id <- rownames(progeny_genos)[i] progeny_vec <- progeny_genos[i, ] - progeny_mat <- base::matrix(progeny_vec, - nrow = base::nrow(parent_genos), - ncol = base::ncol(parent_genos), + progeny_mat <- matrix(progeny_vec, + nrow = nrow(parent_genos), + ncol = ncol(parent_genos), byrow = TRUE) - mismatches <- base::rowSums(parent_genos != progeny_mat, na.rm = TRUE) - comparisons <- base::rowSums(!base::is.na(parent_genos) & !base::is.na(progeny_mat)) + mismatches <- rowSums(parent_genos != progeny_mat, na.rm = TRUE) + comparisons <- rowSums(!is.na(parent_genos) & !is.na(progeny_mat)) percent_mismatch <- (mismatches / comparisons) * 100 - percent_mismatch[base::is.nan(percent_mismatch)] <- NA + percent_mismatch[is.nan(percent_mismatch)] <- NA if (exclude_self_match) { - self_idx <- base::rownames(parent_genos) == prog_id + self_idx <- rownames(parent_genos) == prog_id percent_mismatch[self_idx] <- NA_real_ } if (top_n > 1) { - valid_idx <- base::which(!base::is.na(percent_mismatch)) - ranked_idx <- valid_idx[base::order(percent_mismatch[valid_idx], + valid_idx <- which(!is.na(percent_mismatch)) + ranked_idx <- valid_idx[order(percent_mismatch[valid_idx], -comparisons[valid_idx])] - top_idx <- ranked_idx[base::seq_len(base::min(top_n, base::length(ranked_idx)))] - best_idx <- top_idx[base::seq_len(base::min(1L, base::length(top_idx)))] + top_idx <- ranked_idx[seq_len(min(top_n, length(ranked_idx)))] + best_idx <- top_idx[seq_len(min(1L, length(top_idx)))] } else { - best_idx <- base::which.min(percent_mismatch) + best_idx <- which.min(percent_mismatch) } - if (base::length(best_idx) == 0) { + if (length(best_idx) == 0) { data.table::set(results_dt, i, "markers_tested", 0L) data.table::set(results_dt, i, "status", "low_markers") next } best_markers <- comparisons[best_idx] - best_error <- base::round(percent_mismatch[best_idx], 2) - data.table::set(results_dt, i, "best_match", base::rownames(parent_genos)[best_idx]) + best_error <- round(percent_mismatch[best_idx], 2) + data.table::set(results_dt, i, "best_match", rownames(parent_genos)[best_idx]) data.table::set(results_dt, i, "mendelian_error_pct", best_error) - data.table::set(results_dt, i, "markers_tested", base::as.integer(best_markers)) + data.table::set(results_dt, i, "markers_tested", as.integer(best_markers)) data.table::set(results_dt, i, "status", assign_status(best_markers, best_error)) data.table::set(results_dt, i, "tied_candidates", - base::as.integer(base::sum(percent_mismatch == percent_mismatch[best_idx], + as.integer(sum(percent_mismatch == percent_mismatch[best_idx], na.rm = TRUE))) - if (top_n > 1 && base::length(top_idx) > 1) { - rank_row <- base::list(id = prog_id) - for (k in base::seq(2, base::length(top_idx))) { - rank_row[[base::paste0("best_match_", k)]] <- base::rownames(parent_genos)[top_idx[k]] - rank_row[[base::paste0("mendelian_error_pct_", k)]] <- base::round(percent_mismatch[top_idx[k]], 2) - rank_row[[base::paste0("markers_tested_", k)]] <- base::as.integer(comparisons[top_idx[k]]) + if (top_n > 1 && length(top_idx) > 1) { + rank_row <- list(id = prog_id) + for (k in seq(2, length(top_idx))) { + rank_row[[paste0("best_match_", k)]] <- rownames(parent_genos)[top_idx[k]] + rank_row[[paste0("mendelian_error_pct_", k)]] <- round(percent_mismatch[top_idx[k]], 2) + rank_row[[paste0("markers_tested_", k)]] <- as.integer(comparisons[top_idx[k]]) } extra_list[[i]] <- data.table::as.data.table(rank_row) } } # Add ranks 2..top_n as suffix columns, keeping the progeny order - extra_dt <- data.table::rbindlist(base::Filter(Negate(base::is.null), extra_list), + extra_dt <- data.table::rbindlist(Filter(Negate(is.null), extra_list), fill = TRUE, use.names = TRUE) if (top_n > 1) { - for (k in base::seq(2, top_n)) { - cols <- base::paste0(c("best_match_", "mendelian_error_pct_", "markers_tested_"), k) + for (k in seq(2, top_n)) { + cols <- paste0(c("best_match_", "mendelian_error_pct_", "markers_tested_"), k) for (col in cols) { - vals <- if (base::nrow(extra_dt) > 0 && col %in% base::names(extra_dt)) - extra_dt[[col]][base::match(results_dt$id, extra_dt$id)] else - base::rep(if (grepl("^best_match", col)) NA_character_ else - if (grepl("^markers", col)) NA_integer_ else NA_real_, base::nrow(results_dt)) + vals <- if (nrow(extra_dt) > 0 && col %in% names(extra_dt)) + extra_dt[[col]][match(results_dt$id, extra_dt$id)] else + rep(if (grepl("^best_match", col)) NA_character_ else + if (grepl("^markers", col)) NA_integer_ else NA_real_, nrow(results_dt)) data.table::set(results_dt, j = col, value = vals) } } @@ -360,38 +359,38 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, female_parent = female_parent_candidates$id) if (!allow_parent_selfing) { parent_pairs <- parent_pairs[male_parent != female_parent] - if (verbose) base::cat("Parent selfing is disallowed. Pairs with identical parents are removed.\n") + if (verbose) cat("Parent selfing is disallowed. Pairs with identical parents are removed.\n") } - if (base::nrow(parent_pairs) == 0) stop("No valid parent pairs to test.") + if (nrow(parent_pairs) == 0) stop("No valid parent pairs to test.") male_parent_genos_mat <- genos_mat[parent_pairs$male_parent, , drop = FALSE] female_parent_genos_mat <- genos_mat[parent_pairs$female_parent, , drop = FALSE] progeny_ids <- progeny_candidates$id - n_all <- base::length(progeny_ids) + n_all <- length(progeny_ids) if (pedigree_mode) { - kmale <- base::ifelse(is_unknown(progeny_candidates$male_parent), + kmale <- ifelse(is_unknown(progeny_candidates$male_parent), NA_character_, progeny_candidates$male_parent) - kfem <- base::ifelse(is_unknown(progeny_candidates$female_parent), + kfem <- ifelse(is_unknown(progeny_candidates$female_parent), NA_character_, progeny_candidates$female_parent) } else { - kmale <- kfem <- base::rep(NA_character_, n_all) + kmale <- kfem <- rep(NA_character_, n_all) } - search_mode <- base::ifelse(!base::is.na(kmale), "known_male", - base::ifelse(!base::is.na(kfem), "known_female", "pair")) - pair_idx <- base::which(search_mode == "pair") + search_mode <- ifelse(!is.na(kmale), "known_male", + ifelse(!is.na(kfem), "known_female", "pair")) + pair_idx <- which(search_mode == "pair") # Full pair matrices are only built for progeny with no known parent progeny_mat <- genos_mat[progeny_ids[pair_idx], , drop = FALSE] - n_progeny <- base::nrow(progeny_mat) - n_pairs <- base::nrow(parent_pairs) + n_progeny <- nrow(progeny_mat) + n_pairs <- nrow(parent_pairs) - mismatch_mat <- base::matrix( - base::vapply(base::seq_len(n_progeny), function(j) { + mismatch_mat <- matrix( + vapply(seq_len(n_progeny), function(j) { progeny_vec <- progeny_mat[j, ] - progeny_pair_mat <- base::matrix(progeny_vec, + progeny_pair_mat <- matrix(progeny_vec, nrow = n_pairs, - ncol = base::ncol(progeny_mat), + ncol = ncol(progeny_mat), byrow = TRUE) - base::rowSums( + rowSums( .mend_mismatch(male_parent_genos_mat, female_parent_genos_mat, progeny_pair_mat, ploidy), na.rm = TRUE @@ -400,14 +399,14 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, nrow = n_pairs, ncol = n_progeny ) - comparison_mat <- base::matrix( - base::vapply(base::seq_len(n_progeny), function(j) { + comparison_mat <- matrix( + vapply(seq_len(n_progeny), function(j) { progeny_vec <- progeny_mat[j, ] - progeny_pair_mat <- base::matrix(progeny_vec, + progeny_pair_mat <- matrix(progeny_vec, nrow = n_pairs, - ncol = base::ncol(progeny_mat), + ncol = ncol(progeny_mat), byrow = TRUE) - base::rowSums( + rowSums( .mend_testable(male_parent_genos_mat, female_parent_genos_mat, progeny_pair_mat, ploidy) ) @@ -416,7 +415,7 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, ) pct_mismatch_mat <- (mismatch_mat / comparison_mat) * 100 - pct_mismatch_mat[base::is.nan(pct_mismatch_mat)] <- NA + pct_mismatch_mat[is.nan(pct_mismatch_mat)] <- NA results_dt <- data.table::data.table( id = progeny_ids, @@ -429,19 +428,19 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, ) # Known-parent bookkeeping (pedigree mode only) - known_id_vec <- base::ifelse(search_mode == "known_male", kmale, - base::ifelse(search_mode == "known_female", kfem, + known_id_vec <- ifelse(search_mode == "known_male", kmale, + ifelse(search_mode == "known_female", kfem, NA_character_)) if (pedigree_mode) { - known_err <- base::rep(NA_real_, n_all) - for (j in base::which(search_mode != "pair")) { + known_err <- rep(NA_real_, n_all) + for (j in which(search_mode != "pair")) { kid <- known_id_vec[j] - if (kid %in% base::rownames(genos_hom_mat)) { + if (kid %in% rownames(genos_hom_mat)) { kh <- genos_hom_mat[kid, ] ph <- genos_hom_mat[progeny_ids[j], ] - n_cmp <- base::sum(!base::is.na(kh) & !base::is.na(ph)) + n_cmp <- sum(!is.na(kh) & !is.na(ph)) known_err[j] <- if (n_cmp == 0) NA_real_ else - base::round(base::sum(kh != ph, na.rm = TRUE) / n_cmp * 100, 2) + round(sum(kh != ph, na.rm = TRUE) / n_cmp * 100, 2) } } data.table::set(results_dt, j = "search_mode", value = search_mode) @@ -452,19 +451,19 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, "female_parent", known_id_vec[search_mode == "known_female"]) } - results_list <- base::vector("list", n_all) - for (j in base::seq_len(n_all)) { + results_list <- vector("list", n_all) + for (j in seq_len(n_all)) { prog_id <- progeny_ids[j] mode_j <- search_mode[j] if (mode_j == "pair") { pairs_j <- parent_pairs - pj <- base::match(j, pair_idx) + pj <- match(j, pair_idx) percent_mismatch <- pct_mismatch_mat[, pj] comparisons <- comparison_mat[, pj] } else { known_id <- known_id_vec[j] - if (!known_id %in% base::rownames(genos_mat)) { + if (!known_id %in% rownames(genos_mat)) { warning("Known parent '", known_id, "' of progeny '", prog_id, "' is not in the genotype file; this progeny cannot be tested.", call. = FALSE) @@ -476,8 +475,8 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, male_parent_candidates$id if (!allow_parent_selfing) cand <- cand[cand != known_id] if (exclude_self_match) cand <- cand[cand != prog_id] - cand <- cand[cand %in% base::rownames(genos_mat)] - if (base::length(cand) == 0) { + cand <- cand[cand %in% rownames(genos_mat)] + if (length(cand) == 0) { data.table::set(results_dt, j, "markers_tested", 0L) data.table::set(results_dt, j, "status", "low_markers") next @@ -487,12 +486,12 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, data.table::data.table(male_parent = cand, female_parent = known_id) pm <- genos_mat[pairs_j$male_parent, , drop = FALSE] pf <- genos_mat[pairs_j$female_parent, , drop = FALSE] - pg <- base::matrix(genos_mat[prog_id, ], nrow = base::nrow(pairs_j), - ncol = base::ncol(genos_mat), byrow = TRUE) - mm <- base::rowSums(.mend_mismatch(pm, pf, pg, ploidy), na.rm = TRUE) - comparisons <- base::rowSums(.mend_testable(pm, pf, pg, ploidy)) + pg <- matrix(genos_mat[prog_id, ], nrow = nrow(pairs_j), + ncol = ncol(genos_mat), byrow = TRUE) + mm <- rowSums(.mend_mismatch(pm, pf, pg, ploidy), na.rm = TRUE) + comparisons <- rowSums(.mend_testable(pm, pf, pg, ploidy)) percent_mismatch <- (mm / comparisons) * 100 - percent_mismatch[base::is.nan(percent_mismatch)] <- NA + percent_mismatch[is.nan(percent_mismatch)] <- NA } if (exclude_self_match) { @@ -501,9 +500,9 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, percent_mismatch[self_pair_idx] <- NA_real_ } - min_mismatch_val <- base::min(percent_mismatch, na.rm = TRUE) + min_mismatch_val <- min(percent_mismatch, na.rm = TRUE) - if (base::is.infinite(min_mismatch_val)) { + if (is.infinite(min_mismatch_val)) { data.table::set(results_dt, j, "markers_tested", 0L) data.table::set(results_dt, j, "status", "low_markers") next @@ -511,64 +510,64 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, if (top_n > 1) { # Rank all candidate pairs by error, then by markers tested; keep top_n - valid_idx <- base::which(!base::is.na(percent_mismatch)) - ranked_idx <- valid_idx[base::order(percent_mismatch[valid_idx], + valid_idx <- which(!is.na(percent_mismatch)) + ranked_idx <- valid_idx[order(percent_mismatch[valid_idx], -comparisons[valid_idx])] - best_indices <- ranked_idx[base::seq_len(base::min(top_n, base::length(ranked_idx)))] + best_indices <- ranked_idx[seq_len(min(top_n, length(ranked_idx)))] } else { - best_indices <- base::which(percent_mismatch == min_mismatch_val) - if (base::length(best_indices) > 1) { + best_indices <- which(percent_mismatch == min_mismatch_val) + if (length(best_indices) > 1) { best_markers_per_pair <- comparisons[best_indices] - max_markers <- base::max(best_markers_per_pair) + max_markers <- max(best_markers_per_pair) best_indices <- best_indices[best_markers_per_pair == max_markers] } } best_pairs <- pairs_j[best_indices] best_markers <- comparisons[best_indices[1]] - best_error <- base::round(min_mismatch_val, 2) + best_error <- round(min_mismatch_val, 2) a_status <- assign_status(best_markers, best_error) - if (top_n == 1 && !show_ties && base::nrow(best_pairs) > 1) { - warning("Progeny '", prog_id, "' has ", base::nrow(best_pairs), + if (top_n == 1 && !show_ties && nrow(best_pairs) > 1) { + warning("Progeny '", prog_id, "' has ", nrow(best_pairs), " tied best pairs. Only one is reported as show_ties=FALSE.", call. = FALSE) } report_all <- show_ties || top_n > 1 - num_to_report <- base::min(base::nrow(best_pairs), - if (report_all) base::nrow(best_pairs) else 1) + num_to_report <- min(nrow(best_pairs), + if (report_all) nrow(best_pairs) else 1) data.table::set(results_dt, j, "male_parent", best_pairs$male_parent[1]) data.table::set(results_dt, j, "female_parent", best_pairs$female_parent[1]) - data.table::set(results_dt, j, "mendelian_error_pct", base::sprintf("%.2f", min_mismatch_val)) - data.table::set(results_dt, j, "markers_tested", base::as.integer(best_markers)) + data.table::set(results_dt, j, "mendelian_error_pct", sprintf("%.2f", min_mismatch_val)) + data.table::set(results_dt, j, "markers_tested", as.integer(best_markers)) data.table::set(results_dt, j, "status", a_status) data.table::set(results_dt, j, "tied_candidates", - base::as.integer(base::sum(percent_mismatch == min_mismatch_val, + as.integer(sum(percent_mismatch == min_mismatch_val, na.rm = TRUE))) if (report_all && num_to_report > 1) { - tie_row <- base::list(id = prog_id) - for (k in base::seq(2, num_to_report)) { - tie_row[[base::paste0("male_parent_", k)]] <- best_pairs$male_parent[k] - tie_row[[base::paste0("female_parent_", k)]] <- best_pairs$female_parent[k] - tie_row[[base::paste0("mendelian_error_pct_", k)]] <- - if (top_n > 1) base::round(percent_mismatch[best_indices[k]], 2) else min_mismatch_val - tie_row[[base::paste0("markers_tested_", k)]] <- comparisons[best_indices[k]] + tie_row <- list(id = prog_id) + for (k in seq(2, num_to_report)) { + tie_row[[paste0("male_parent_", k)]] <- best_pairs$male_parent[k] + tie_row[[paste0("female_parent_", k)]] <- best_pairs$female_parent[k] + tie_row[[paste0("mendelian_error_pct_", k)]] <- + if (top_n > 1) round(percent_mismatch[best_indices[k]], 2) else min_mismatch_val + tie_row[[paste0("markers_tested_", k)]] <- comparisons[best_indices[k]] } results_list[[j]] <- data.table::as.data.table(tie_row) } } tie_rows <- data.table::rbindlist( - base::Filter(Negate(base::is.null), results_list), + Filter(Negate(is.null), results_list), fill = TRUE, use.names = TRUE ) - if (base::nrow(tie_rows) > 0) { + if (nrow(tie_rows) > 0) { final_df <- merge(results_dt, tie_rows, by = "id", all.x = TRUE) - for (col in base::names(final_df)) + for (col in names(final_df)) data.table::set(final_df, which(final_df[[col]] == ""), col, NA_character_) } else { final_df <- results_dt @@ -586,40 +585,39 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, #### Verbose output #### if (verbose) { - total_progeny <- base::nrow(final_df) - base::cat("\n=== Parentage Assignment Report ===\n") - base::cat("\nTotal progeny evaluated:", total_progeny, "\n") - base::cat("Method:", method_reported, " | ", - "Error threshold:", error_threshold, "% | ", - "Minimum markers:", min_markers, "\n") + total_progeny <- nrow(final_df) + cat("\n=== Parentage Assignment Report ===\n") + cat(sprintf("\nTotal progeny evaluated: %d\n", total_progeny)) + cat(sprintf("Method: %s | Error threshold: %s%% | Minimum markers: %s\n", + method_reported, error_threshold, min_markers)) - n_pass <- base::nrow(output_list$pass) + n_pass <- nrow(output_list$pass) if (n_pass > 0) { - base::cat(base::sprintf("\n%d progeny passed (%.1f%%).\n", + cat(sprintf("\n%d progeny passed (%.1f%%).\n", n_pass, (n_pass / total_progeny) * 100)) } else { - base::cat("\nNo progeny passed.\n") + cat("\nNo progeny passed.\n") } - n_high <- base::nrow(output_list$high_error) + n_high <- nrow(output_list$high_error) if (n_high > 0) { - base::cat(base::sprintf("\n%d progeny flagged high_error (%.1f%%):\n", + cat(sprintf("\n%d progeny flagged high_error (%.1f%%):\n", n_high, (n_high / total_progeny) * 100)) - base::print(output_list$high_error) + print(output_list$high_error) } else { - base::cat("\nNo progeny flagged for high error.\n") + cat("\nNo progeny flagged for high error.\n") } - n_low <- base::nrow(output_list$low_markers) + n_low <- nrow(output_list$low_markers) if (n_low > 0) { - base::cat(base::sprintf("\n%d progeny flagged low_markers (%.1f%%):\n", + cat(sprintf("\n%d progeny flagged low_markers (%.1f%%):\n", n_low, (n_low / total_progeny) * 100)) - base::print(output_list$low_markers) + print(output_list$low_markers) } else { - base::cat("\nNo progeny flagged for low marker count.\n") + cat("\nNo progeny flagged for low marker count.\n") } - base::cat("\nFull results are included in the returned list as $full_results.\n") + cat("\nFull results are included in the returned list as $full_results.\n") } #### Plot Results #### @@ -629,20 +627,20 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, call. = FALSE) } else { plot_df <- final_df[!is.na(final_df$mendelian_error_pct)] - plot_df$mendelian_error_pct <- base::as.numeric(plot_df$mendelian_error_pct) - plot_df$plot_status <- base::ifelse( + plot_df$mendelian_error_pct <- as.numeric(plot_df$mendelian_error_pct) + plot_df$plot_status <- ifelse( plot_df$status == "pass", "pass", - base::ifelse( + ifelse( plot_df$status == "high_error", "high_error", - base::ifelse( + ifelse( plot_df$status == "low_markers", "low_markers", "other"))) - n_total <- base::nrow(plot_df) - n_pass <- base::sum(plot_df$status == "pass", na.rm = TRUE) - n_high <- base::sum(plot_df$status == "high_error", na.rm = TRUE) - n_low <- base::sum(plot_df$status == "low_markers", na.rm = TRUE) + n_total <- nrow(plot_df) + n_pass <- sum(plot_df$status == "pass", na.rm = TRUE) + n_high <- sum(plot_df$status == "high_error", na.rm = TRUE) + n_low <- sum(plot_df$status == "low_markers", na.rm = TRUE) - threshold_label <- base::paste0( + threshold_label <- paste0( "Error Threshold: ", error_threshold, "% | ", "Pass: ", n_pass, " | ", "High Error: ", n_high, " | ", @@ -670,7 +668,7 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, ) + ggplot2::labs( title = "Parentage Mendelian Error Distribution", - subtitle = base::paste0("Progeny Tested: ", n_total, + subtitle = paste0("Progeny Tested: ", n_total, "\n \n", threshold_label), x = "Mendelian Error (%)", y = "Number of Progeny", @@ -679,10 +677,10 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, ggplot2::theme_classic(base_size = 13) + ggplot2::theme(legend.position = "top") - base::print(p) + print(p) output_list$plot <- p } } - return(base::invisible(output_list)) + return(invisible(output_list)) } diff --git a/R/ped_to_dosage.R b/R/ped_to_dosage.R index 14faa92..232ebd7 100644 --- a/R/ped_to_dosage.R +++ b/R/ped_to_dosage.R @@ -5,7 +5,7 @@ #' columns, coded as the number of copies of the counted allele (0, 1, 2) #' with NA for missing calls. PLINK `.ped` files are diploid only. #' -#' The `.ped` format stores allele letters rather than a reference/alternate +#' The `.ped` format stores allele letters and has no reference/alternate #' pair, so the counted allele (allele B) is chosen per marker. By default it #' is the alphabetically / numerically last allele observed at the marker, so #' `1/2` coded files count allele `2` and an `A/G` marker counts `G`. To code a @@ -71,112 +71,112 @@ ped_to_dosage <- function(ped, verbose = TRUE) { #### Input Validation #### - format <- base::match.arg(format) - if (!base::is.character(ped) || base::length(ped) != 1) + format <- match.arg(format) + if (!is.character(ped) || length(ped) != 1) stop("ped must be a path to a PLINK .ped file.") - if (!base::file.exists(ped)) + if (!file.exists(ped)) stop("PED file not found: ", ped) - if (!base::is.null(counted_allele) && - (!base::is.character(counted_allele) || base::is.null(base::names(counted_allele)))) + if (!is.null(counted_allele) && + (!is.character(counted_allele) || is.null(names(counted_allele)))) stop("counted_allele must be a named character vector (names = marker names).") #### Read the .ped #### raw <- data.table::fread(ped, header = FALSE, colClasses = "character", data.table = FALSE) - n_allele_cols <- base::ncol(raw) - 6L + n_allele_cols <- ncol(raw) - 6L if (n_allele_cols < 2 || n_allele_cols %% 2 != 0) stop("The .ped file must have 6 pedigree columns followed by two allele ", "columns per marker (compound genotypes are not supported).") n_markers <- n_allele_cols %/% 2L ids <- raw[[2]] - dup <- base::unique(ids[base::duplicated(ids)]) - if (base::length(dup) > 0) + dup <- unique(ids[duplicated(ids)]) + if (length(dup) > 0) stop("Duplicated individual IDs (column 2) in the .ped file: ", - base::paste(dup, collapse = ", ")) + paste(dup, collapse = ", ")) - alleles <- base::as.matrix(raw[, -(1:6), drop = FALSE]) + alleles <- as.matrix(raw[, -(1:6), drop = FALSE]) alleles[alleles %in% missing_codes] <- NA - a1 <- alleles[, base::seq(1, n_allele_cols, by = 2), drop = FALSE] - a2 <- alleles[, base::seq(2, n_allele_cols, by = 2), drop = FALSE] - miss <- base::is.na(a1) | base::is.na(a2) + a1 <- alleles[, seq(1, n_allele_cols, by = 2), drop = FALSE] + a2 <- alleles[, seq(2, n_allele_cols, by = 2), drop = FALSE] + miss <- is.na(a1) | is.na(a2) a1[miss] <- NA a2[miss] <- NA #### Marker names #### - map_path <- if (base::identical(map, "auto")) { - cand <- base::sub("\\.ped$", ".map", ped, ignore.case = TRUE) - if (cand != ped && base::file.exists(cand)) cand else NULL + map_path <- if (identical(map, "auto")) { + cand <- sub("\\.ped$", ".map", ped, ignore.case = TRUE) + if (cand != ped && file.exists(cand)) cand else NULL } else { map } - if (!base::is.null(map_path)) { - if (!base::file.exists(map_path)) + if (!is.null(map_path)) { + if (!file.exists(map_path)) stop("MAP file not found: ", map_path) map_df <- data.table::fread(map_path, header = FALSE, colClasses = "character", data.table = FALSE) - if (base::nrow(map_df) != n_markers) - stop("The .map file lists ", base::nrow(map_df), " markers but the .ped file has ", + if (nrow(map_df) != n_markers) + stop("The .map file lists ", nrow(map_df), " markers but the .ped file has ", n_markers, ".") snp_names <- map_df[[2]] - no_name <- base::is.na(snp_names) | snp_names %in% c("", ".") - snp_names[no_name] <- base::paste(map_df[[1]][no_name], - map_df[[base::ncol(map_df)]][no_name], sep = "_") + no_name <- is.na(snp_names) | snp_names %in% c("", ".") + snp_names[no_name] <- paste(map_df[[1]][no_name], + map_df[[ncol(map_df)]][no_name], sep = "_") } else { - snp_names <- base::paste0("SNP", base::seq_len(n_markers)) + snp_names <- paste0("SNP", seq_len(n_markers)) } - if (base::any(base::duplicated(snp_names))) { + if (any(duplicated(snp_names))) { warning("Duplicated marker names were made unique with a numeric suffix.", call. = FALSE) - snp_names <- base::make.unique(snp_names, sep = "_") + snp_names <- make.unique(snp_names, sep = "_") } #### Counted allele per marker #### - both <- base::rbind(a1, a2) - n_alleles <- base::apply(both, 2, function(v) base::length(base::unique(v[!base::is.na(v)]))) - counted <- base::apply(both, 2, function(v) { - v <- v[!base::is.na(v)] - if (base::length(v) == 0) NA_character_ else base::max(v) + both <- rbind(a1, a2) + n_alleles <- apply(both, 2, function(v) length(unique(v[!is.na(v)]))) + counted <- apply(both, 2, function(v) { + v <- v[!is.na(v)] + if (length(v) == 0) NA_character_ else max(v) }) - base::names(counted) <- snp_names + names(counted) <- snp_names - if (!base::is.null(counted_allele)) { - use <- snp_names[snp_names %in% base::names(counted_allele)] + if (!is.null(counted_allele)) { + use <- snp_names[snp_names %in% names(counted_allele)] counted[use] <- counted_allele[use] } #### Drop multiallelic markers #### keep <- n_alleles <= 2 - if (base::any(!keep)) - warning(base::sum(!keep), " marker(s) with more than two alleles were dropped.", + if (any(!keep)) + warning(sum(!keep), " marker(s) with more than two alleles were dropped.", call. = FALSE) - if (base::sum(keep) == 0) + if (sum(keep) == 0) stop("No biallelic markers remain in the .ped file.") #### Dosage = copies of the counted allele #### - cnt <- base::matrix(counted, nrow = base::nrow(a1), ncol = n_markers, byrow = TRUE) + cnt <- matrix(counted, nrow = nrow(a1), ncol = n_markers, byrow = TRUE) geno <- (a1 == cnt) + (a2 == cnt) - base::storage.mode(geno) <- "integer" - base::dimnames(geno) <- base::list(ids, snp_names) + storage.mode(geno) <- "integer" + dimnames(geno) <- list(ids, snp_names) geno <- geno[, keep, drop = FALSE] counted <- counted[keep] #### Verbose report #### if (verbose) { - base::cat("\n=== PED Conversion Report ===\n") - base::cat("Ploidy: 2 (PLINK .ped)\n") - base::cat("Marker names from: ", if (base::is.null(map_path)) "position (no .map)" else map_path, "\n") - base::cat("Markers read: ", n_markers, "\n") - base::cat("Markers returned: ", base::ncol(geno), "\n") - base::cat("Individuals returned: ", base::nrow(geno), "\n") - base::cat(base::sprintf("Missing calls: %.2f%%\n", - 100 * base::mean(base::is.na(geno)))) + cat("\n=== PED Conversion Report ===\n") + cat("Ploidy: 2 (PLINK .ped)\n") + cat(sprintf("Marker names from: %s\n", + if (is.null(map_path)) "position (no .map)" else map_path)) + cat(sprintf("Markers read: %s\n", n_markers)) + cat(sprintf("Markers returned: %d\n", ncol(geno))) + cat(sprintf("Individuals returned: %d\n", nrow(geno))) + cat(sprintf("Missing calls: %.2f%%\n", 100 * mean(is.na(geno)))) } #### Return in the requested layout #### if (format == "matrix") { - base::attr(geno, "ploidy") <- 2L - base::attr(geno, "counted_allele") <- counted + attr(geno, "ploidy") <- 2L + attr(geno, "counted_allele") <- counted return(geno) } diff --git a/R/utils.R b/R/utils.R index c5efe63..ca73095 100644 --- a/R/utils.R +++ b/R/utils.R @@ -25,8 +25,8 @@ utils::globalVariables(c( #' @param ploidy even integer ploidy level. #' @return Object matching \code{g} holding the gamete dosage bound. #' @noRd -.gamete_lo <- function(g, ploidy) base::pmax(0L, g - ploidy / 2L) -.gamete_hi <- function(g, ploidy) base::pmin(ploidy / 2L, g) +.gamete_lo <- function(g, ploidy) pmax(0L, g - ploidy / 2L) +.gamete_hi <- function(g, ploidy) pmin(ploidy / 2L, g) #' Flag impossible offspring dosages #' @@ -79,10 +79,10 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @noRd .mend_testable <- function(male, female, offspring, ploidy) { if (ploidy %% 2 == 0) - return(!base::is.na(male) & !base::is.na(female) & !base::is.na(offspring)) - o_hom <- !base::is.na(offspring) & (offspring == 0 | offspring == ploidy) - m_hom <- !base::is.na(male) & (male == 0 | male == ploidy) - f_hom <- !base::is.na(female) & (female == 0 | female == ploidy) + return(!is.na(male) & !is.na(female) & !is.na(offspring)) + o_hom <- !is.na(offspring) & (offspring == 0 | offspring == ploidy) + m_hom <- !is.na(male) & (male == 0 | male == ploidy) + f_hom <- !is.na(female) & (female == 0 | female == ploidy) o_hom & (m_hom | f_hom) } @@ -94,10 +94,10 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @return Invisibly TRUE if valid; otherwise an error is thrown. #' @noRd .check_ploidy <- function(ploidy) { - if (base::length(ploidy) != 1 || !base::is.numeric(ploidy) || - base::is.na(ploidy) || ploidy < 2 || ploidy != base::round(ploidy)) - base::stop("ploidy must be an integer >= 2.") - base::invisible(TRUE) + if (length(ploidy) != 1 || !is.numeric(ploidy) || + is.na(ploidy) || ploidy < 2 || ploidy != round(ploidy)) + stop("ploidy must be an integer >= 2.") + invisible(TRUE) } #### Column name helpers #### @@ -112,7 +112,7 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @return Normalized character vector. #' @noRd .normalize_name <- function(x) { - base::gsub("[ .-]+", "_", base::tolower(base::trimws(x))) + gsub("[ .-]+", "_", tolower(trimws(x))) } #' Rename columns to their standard names, ignoring case @@ -127,23 +127,23 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @return \code{x} with matching columns renamed. #' @noRd .standardize_names <- function(x, cols) { - nm <- base::names(x) + nm <- names(x) norm <- .normalize_name(nm) - old <- new <- base::character(0) + old <- new <- character(0) for (col in cols) { if (col %in% nm) next - hit <- base::which(norm == col & !(nm %in% old)) - if (base::length(hit) >= 1) { + hit <- which(norm == col & !(nm %in% old)) + if (length(hit) >= 1) { old <- c(old, nm[hit[1]]) new <- c(new, col) } } - if (base::length(old) == 0) return(x) + if (length(old) == 0) return(x) if (data.table::is.data.table(x)) { x <- data.table::copy(x) data.table::setnames(x, old, new) } else { - base::names(x)[base::match(old, base::names(x))] <- new + names(x)[match(old, names(x))] <- new } x } @@ -160,15 +160,15 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @return One of "vcf", "plink", "text", "data.frame", or "matrix". #' @noRd .genotype_source <- function(x) { - if (base::inherits(x, "vcfR")) return("vcf") - if (base::is.matrix(x)) return("matrix") - if (base::is.data.frame(x)) return("data.frame") - if (base::is.character(x) && base::length(x) == 1) { - if (base::grepl("\\.vcf(\\.gz)?$", x, ignore.case = TRUE)) return("vcf") - if (base::grepl("\\.ped$", x, ignore.case = TRUE)) return("plink") + if (inherits(x, "vcfR")) return("vcf") + if (is.matrix(x)) return("matrix") + if (is.data.frame(x)) return("data.frame") + if (is.character(x) && length(x) == 1) { + if (grepl("\\.vcf(\\.gz)?$", x, ignore.case = TRUE)) return("vcf") + if (grepl("\\.ped$", x, ignore.case = TRUE)) return("plink") return("text") } - base::stop("Genotypes must be a file path (TXT/TSV/CSV, VCF or PLINK .ped), ", + stop("Genotypes must be a file path (TXT/TSV/CSV, VCF or PLINK .ped), ", "a vcfR object, a data.frame / data.table, or a matrix.") } @@ -196,7 +196,7 @@ mendelian_error <- function(male, female, offspring, ploidy) { format = c("data.frame", "matrix"), verbose = FALSE, counted_allele = NULL) { - format <- base::match.arg(format) + format <- match.arg(format) src <- .genotype_source(x) # VCF: convert straight to the requested layout @@ -205,8 +205,8 @@ mendelian_error <- function(male, female, offspring, ploidy) { # PLINK .ped: diploid only; a .map next to the file supplies marker names if (src == "plink") { - if (!base::is.null(ploidy) && ploidy != 2) - base::stop("PLINK .ped files are diploid; use ploidy = 2.") + if (!is.null(ploidy) && ploidy != 2) + stop("PLINK .ped files are diploid; use ploidy = 2.") return(ped_to_dosage(x, format = format, counted_allele = counted_allele, verbose = verbose)) } @@ -214,39 +214,39 @@ mendelian_error <- function(male, female, offspring, ploidy) { # Text file: read to a table, then fall through to the in-memory handling from_text <- src == "text" if (from_text) { - if (!base::file.exists(x)) - base::stop("Genotype file not found: ", x) + if (!file.exists(x)) + stop("Genotype file not found: ", x) x <- data.table::fread(x, sep = "auto", data.table = FALSE) src <- "data.frame" } # ID column: "id" in any case (id, ID, Id, ...); exact "id" preferred id_col <- if (src == "data.frame") { - nm <- base::names(x) + nm <- names(x) hit <- if ("id" %in% nm) "id" else nm[.normalize_name(nm) == "id"] - if (base::length(hit) >= 1) hit[1] else NA_character_ + if (length(hit) >= 1) hit[1] else NA_character_ } else { NA_character_ } # Text files have no row names, so the first column holds the IDs - if (from_text && format == "matrix" && base::is.na(id_col)) - id_col <- base::names(x)[1] + if (from_text && format == "matrix" && is.na(id_col)) + id_col <- names(x)[1] if (format == "data.frame") { if (src == "matrix") return(data.table::as.data.table(x, keep.rownames = "id")) out <- data.table::as.data.table(x) - if (!base::is.na(id_col) && id_col != "id") + if (!is.na(id_col) && id_col != "id") data.table::setnames(out, id_col, "id") return(out) } # format == "matrix": individuals in named rows, markers in columns - if (src == "matrix" || base::is.na(id_col)) return(x) - x <- base::as.data.frame(x) # data.table subsetting differs from data.frame - ids <- base::as.character(x[[id_col]]) - x <- base::as.matrix(x[, base::names(x) != id_col, drop = FALSE]) - base::rownames(x) <- ids + if (src == "matrix" || is.na(id_col)) return(x) + x <- as.data.frame(x) # data.table subsetting differs from data.frame + ids <- as.character(x[[id_col]]) + x <- as.matrix(x[, names(x) != id_col, drop = FALSE]) + rownames(x) <- ids x } @@ -323,7 +323,7 @@ QPsolve <- function(Y, X) { if (ploidy %% 2 == 0) { h <- ploidy / 2 parent_bad <- function(g) { - !(offspring >= base::pmax(0, g - h) & offspring <= base::pmin(g, h) + h) + !(offspring >= pmax(0, g - h) & offspring <= pmin(g, h) + h) } male_bad <- parent_bad(male) female_bad <- parent_bad(female) @@ -334,11 +334,11 @@ QPsolve <- function(Y, X) { } # missing calls are named after every individual without a call, e.g. # missing_male, missing_male_female, missing_progeny_female - na_p <- base::is.na(offspring); na_m <- base::is.na(male); na_f <- base::is.na(female) - missing_label <- base::paste0("missing", - base::ifelse(na_p, "_progeny", ""), - base::ifelse(na_m, "_male", ""), - base::ifelse(na_f, "_female", "")) + na_p <- is.na(offspring); na_m <- is.na(male); na_f <- is.na(female) + missing_label <- paste0("missing", + ifelse(na_p, "_progeny", ""), + ifelse(na_m, "_male", ""), + ifelse(na_f, "_female", "")) dplyr::case_when( na_p | na_m | na_f ~ missing_label, !tt ~ "untestable", @@ -364,6 +364,10 @@ QPsolve <- function(Y, X) { ) return(unique(out[[1]])) } + # A single string that looks like a path but does not exist is a typo, not an id + if (is.character(x) && length(x) == 1 && + grepl("[/\\\\]|\\.(txt|csv|tsv|tab)$", x, ignore.case = TRUE)) + stop("Could not read founders list. Ensure it is a plain text or CSV/TSV file.") if (is.atomic(x)) return(unique(as.character(x))) stop("founders_file must be a file path, a vector of ids, or a data.frame.") } diff --git a/R/validate_pedigree.R b/R/validate_pedigree.R index 00ce2f5..56805ad 100644 --- a/R/validate_pedigree.R +++ b/R/validate_pedigree.R @@ -178,11 +178,11 @@ validate_pedigree <- function(pedigree_file, genotypes_file, # Column names are matched ignoring case (e.g. ID, Male_Parent, FEMALE_PARENT) pedigree <- .standardize_names(pedigree, c("id", "male_parent", "female_parent")) required_ped_cols <- c("id", "male_parent", "female_parent") - missing_cols <- base::setdiff(required_ped_cols, base::names(pedigree)) - if (base::length(missing_cols) > 0) + missing_cols <- setdiff(required_ped_cols, names(pedigree)) + if (length(missing_cols) > 0) stop("Pedigree file missing required columns: ", - base::paste(missing_cols, collapse = ", ")) - if (!"id" %in% base::names(genos)) + paste(missing_cols, collapse = ", ")) + if (!"id" %in% names(genos)) stop("Genotypes file must have an 'id' column") pedigree[, male_parent := as.character(male_parent)] @@ -193,12 +193,12 @@ validate_pedigree <- function(pedigree_file, genotypes_file, founder_ids <- .read_founders(founders_file) #### Build genotype matrices #### - genos_mat <- base::as.matrix(genos, rownames = "id") + genos_mat <- as.matrix(genos, rownames = "id") genos_hom <- data.table::copy(genos) - marker_cols <- base::setdiff(base::names(genos_hom), "id") + marker_cols <- setdiff(names(genos_hom), "id") for (col in marker_cols) - genos_hom[base::get(col) > 0 & base::get(col) < ploidy, (col) := NA_integer_] - genos_hom_mat <- base::as.matrix(genos_hom, rownames = "id") + genos_hom[get(col) > 0 & get(col) < ploidy, (col) := NA_integer_] + genos_hom_mat <- as.matrix(genos_hom, rownames = "id") #### Identify trios missing from the genotype file #### valid_ids <- as.character(genos$id) @@ -208,20 +208,20 @@ validate_pedigree <- function(pedigree_file, genotypes_file, no_geno_rows <- pedigree[!(id %in% valid_ids) | (!(male_parent %in% valid_ids) & male_parent != "0") | (!(female_parent %in% valid_ids) & female_parent != "0")] - if (base::nrow(no_geno_rows) > 0 && verbose) - base::cat("Found", base::nrow(no_geno_rows), - "trios with missing genotype data; flagged as no_genotype_data.\n") + if (nrow(no_geno_rows) > 0 && verbose) + cat(sprintf("Found %d trios with missing genotype data; flagged as no_genotype_data.\n", + nrow(no_geno_rows))) pedigree <- has_geno - if (base::nrow(pedigree) == 0) + if (nrow(pedigree) == 0) stop("No valid trios remain after filtering for genotype availability.") #### Per-marker mismatch counters (filled in the trio loop) #### - marker_names <- base::colnames(genos_mat) - marker_tested_n <- base::numeric(base::length(marker_names)) - marker_mismatch_n <- base::numeric(base::length(marker_names)) + marker_names <- colnames(genos_mat) + marker_tested_n <- numeric(length(marker_names)) + marker_mismatch_n <- numeric(length(marker_names)) #### Main trio evaluation loop #### - results_list <- base::lapply(base::seq_len(base::nrow(pedigree)), function(i) { + results_list <- lapply(seq_len(nrow(pedigree)), function(i) { prog_id <- pedigree$id[i] male_parent_id <- pedigree$male_parent[i] female_parent_id <- pedigree$female_parent[i] @@ -254,11 +254,11 @@ validate_pedigree <- function(pedigree_file, genotypes_file, progeny_vec, ploidy) tt_vec <- .mend_testable(male_parent_vec, female_parent_vec, progeny_vec, ploidy) - mismatches <- base::sum(mm_vec, na.rm = TRUE) - markers_tested <- base::sum(tt_vec) - marker_tested_n <<- marker_tested_n + base::as.numeric(tt_vec) + mismatches <- sum(mm_vec, na.rm = TRUE) + markers_tested <- sum(tt_vec) + marker_tested_n <<- marker_tested_n + as.numeric(tt_vec) marker_mismatch_n <<- marker_mismatch_n + - base::as.numeric((mm_vec %in% TRUE) & tt_vec) + as.numeric((mm_vec %in% TRUE) & tt_vec) if (markers_tested == 0) { status <- "no_data" correction_decision <- "none" @@ -276,15 +276,15 @@ validate_pedigree <- function(pedigree_file, genotypes_file, progeny_hom <- genos_hom_mat[prog_id, ] male_parent_hom <- genos_hom_mat[male_parent_id, ] female_parent_hom <- genos_hom_mat[female_parent_id, ] - male_comparisons <- base::sum(!base::is.na(male_parent_hom) & - !base::is.na(progeny_hom)) + male_comparisons <- sum(!is.na(male_parent_hom) & + !is.na(progeny_hom)) male_parent_error_pct <- if (male_comparisons == 0) NA_real_ else - base::round((base::sum(male_parent_hom != progeny_hom, na.rm = TRUE) / + round((sum(male_parent_hom != progeny_hom, na.rm = TRUE) / male_comparisons) * 100, 2) - female_comparisons <- base::sum(!base::is.na(female_parent_hom) & - !base::is.na(progeny_hom)) + female_comparisons <- sum(!is.na(female_parent_hom) & + !is.na(progeny_hom)) female_parent_error_pct <- if (female_comparisons == 0) NA_real_ else - base::round((base::sum(female_parent_hom != progeny_hom, na.rm = TRUE) / + round((sum(female_parent_hom != progeny_hom, na.rm = TRUE) / female_comparisons) * 100, 2) male_acceptable <- !is.na(male_parent_error_pct) && male_parent_error_pct <= single_parent_error_threshold @@ -294,7 +294,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, # Neither parent is excluded by the homozygous check. If the trio # still fails (e.g. a recorded self whose offspring carries alleles # absent from the parent), the pedigree is wrong but the faulty - # parent cannot be identified: flag it instead of keeping it. + # parent cannot be identified, so the trio is flagged and the parents are not kept. correction_decision <- if (error_pct > trio_error_threshold) "unresolved" else "keep_both" } else if (male_acceptable && !female_acceptable) { correction_decision <- "remove_female_parent" @@ -314,7 +314,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, id = prog_id, orig_male_parent = male_parent_id, orig_female_parent = female_parent_id, - trio_mendelian_error_pct = base::round(error_pct, 2), + trio_mendelian_error_pct = round(error_pct, 2), trio_markers_tested = markers_tested, status = status, recommended_correction = correction_decision, @@ -326,7 +326,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, final_df <- data.table::rbindlist(results_list) #### Append no_genotype_data rows #### - if (base::nrow(no_geno_rows) > 0) { + if (nrow(no_geno_rows) > 0) { no_geno_df <- data.table::data.table( id = no_geno_rows$id, orig_male_parent = no_geno_rows$male_parent, @@ -343,12 +343,12 @@ validate_pedigree <- function(pedigree_file, genotypes_file, #### Build corrected pedigree #### corrected_pedigree <- data.table::copy(original_pedigree) - for (i in base::seq_len(base::nrow(final_df))) { + for (i in seq_len(nrow(final_df))) { prog_id <- final_df$id[i] decision <- final_df$recommended_correction[i] - row_idx <- base::which(corrected_pedigree$id == prog_id) + row_idx <- which(corrected_pedigree$id == prog_id) # low_markers_* decisions apply the same removal as the plain decision - base_decision <- base::sub("^low_markers_", "", decision) + base_decision <- sub("^low_markers_", "", decision) if (base_decision %in% c("remove_male_parent", "remove_both")) data.table::set(corrected_pedigree, row_idx, "male_parent", "0") if (base_decision %in% c("remove_female_parent", "remove_both")) @@ -357,24 +357,23 @@ validate_pedigree <- function(pedigree_file, genotypes_file, #### Summary output #### if (verbose) { - total_trios <- base::nrow(final_df) - status_counts <- base::table(final_df$status) - base::cat("\n--- Trio Validation Summary ---\n") - base::cat("Total trios in pedigree:", total_trios, "\n") - for (s in base::names(status_counts)) - base::cat(base::sprintf("%-28s: %d (%.1f%%)\n", s, - status_counts[s], - (status_counts[s] / total_trios) * 100)) - base::cat("Error threshold:", trio_error_threshold, "%\n") - base::cat("Homozygous threshold:", single_parent_error_threshold, "%\n") - base::cat("Minimum markers required:", min_markers, "\n\n") - corrections <- base::table(final_df$recommended_correction) - base::cat("Correction summary:\n") - for (decision in base::names(corrections)) + total_trios <- nrow(final_df) + status_counts <- table(final_df$status) + cat("\n--- Trio Validation Summary ---\n") + cat(sprintf("Total trios in pedigree: %d\n", total_trios)) + for (s in names(status_counts)) + cat(sprintf("%-28s: %d (%.1f%%)\n", s, status_counts[s], + (status_counts[s] / total_trios) * 100)) + cat(sprintf("Error threshold: %s%%\n", trio_error_threshold)) + cat(sprintf("Homozygous threshold: %s%%\n", single_parent_error_threshold)) + cat(sprintf("Minimum markers required: %s\n\n", min_markers)) + corrections <- table(final_df$recommended_correction) + cat("Correction summary:\n") + for (decision in names(corrections)) if (decision != "none") - base::cat(" ", decision, ":", corrections[decision], "\n") - base::cat("\n") - base::print(final_df) + cat(sprintf(" %s: %d\n", decision, corrections[decision])) + cat("\n") + print(final_df) } #### Plot results #### @@ -383,20 +382,21 @@ validate_pedigree <- function(pedigree_file, genotypes_file, if (!requireNamespace("ggplot2", quietly = TRUE)) { warning("ggplot2 is required for plot_results = TRUE. Please install it.", call. = FALSE) } else { - plot_df <- final_df[!is.na(final_df$trio_mendelian_error_pct)] # Colour by the trio error first, so no trio above the threshold is shown as pass - plot_df$plot_status <- dplyr::case_when( - plot_df$trio_mendelian_error_pct <= trio_error_threshold ~ "pass", - plot_df$recommended_correction %in% c("remove_male_parent", - "remove_female_parent", - "low_markers_remove_male_parent", - "low_markers_remove_female_parent") ~ "fail_one_parent", - plot_df$recommended_correction %in% c("remove_both", - "low_markers_remove_both") ~ "fail_both_parents", - plot_df$recommended_correction %in% c("unresolved", - "low_markers_unresolved") ~ "fail_unresolved", - TRUE ~ "other" - ) + plot_df <- final_df %>% + dplyr::filter(!is.na(trio_mendelian_error_pct)) %>% + dplyr::mutate(plot_status = dplyr::case_when( + trio_mendelian_error_pct <= trio_error_threshold ~ "pass", + recommended_correction %in% c("remove_male_parent", + "remove_female_parent", + "low_markers_remove_male_parent", + "low_markers_remove_female_parent") ~ "fail_one_parent", + recommended_correction %in% c("remove_both", + "low_markers_remove_both") ~ "fail_both_parents", + recommended_correction %in% c("unresolved", + "low_markers_unresolved") ~ "fail_unresolved", + TRUE ~ "other" + )) n_total <- nrow(plot_df) n_fail <- sum(plot_df$trio_mendelian_error_pct > trio_error_threshold) n_pass <- sum(plot_df$trio_mendelian_error_pct <= trio_error_threshold) @@ -444,29 +444,30 @@ validate_pedigree <- function(pedigree_file, genotypes_file, marker_trio_out <- NULL if (marker_trio_table) { full_trios <- pedigree[male_parent != "0" & female_parent != "0"] - status_list <- base::lapply(base::seq_len(base::nrow(full_trios)), function(i) + status_list <- lapply(seq_len(nrow(full_trios)), function(i) .marker_trio_status(genos_mat[full_trios$male_parent[i], ], genos_mat[full_trios$female_parent[i], ], genos_mat[full_trios$id[i], ], ploidy)) - base::names(status_list) <- base::make.unique( - base::paste(full_trios$id, full_trios$male_parent, full_trios$female_parent, sep = "_")) + names(status_list) <- make.unique( + paste(full_trios$id, full_trios$male_parent, full_trios$female_parent, sep = "_")) marker_trio_out <- data.table::as.data.table( - base::c(base::list(marker = marker_names), status_list)) + c(list(marker = marker_names), status_list)) } #### Per-marker mismatch summary #### marker_summary <- data.table::data.table( marker = marker_names, - trios_tested = base::as.integer(marker_tested_n), - trios_mismatch = base::as.integer(marker_mismatch_n) - ) - marker_summary[, mismatch_pct := base::ifelse( - trios_tested > 0, base::round(trios_mismatch / trios_tested * 100, 2), NA_real_)] - marker_summary <- marker_summary[base::order(-mismatch_pct, -trios_mismatch, marker)] + trios_tested = as.integer(marker_tested_n), + trios_mismatch = as.integer(marker_mismatch_n) + ) %>% + dplyr::mutate(mismatch_pct = ifelse(trios_tested > 0, + round(trios_mismatch / trios_tested * 100, 2), + NA_real_)) %>% + dplyr::arrange(dplyr::desc(mismatch_pct), dplyr::desc(trios_mismatch), marker) #### Compile and return named list #### - output_list <- base::list( + output_list <- list( pass = final_df[status == "pass"], fail = final_df[status == "fail"], low_markers = final_df[status == "low_markers"], @@ -481,5 +482,5 @@ validate_pedigree <- function(pedigree_file, genotypes_file, marker_trio_table = marker_trio_out, plot = p ) - return(base::invisible(output_list)) + return(invisible(output_list)) } diff --git a/R/vcf_to_dosage.R b/R/vcf_to_dosage.R index 7d753cf..149f945 100644 --- a/R/vcf_to_dosage.R +++ b/R/vcf_to_dosage.R @@ -107,27 +107,27 @@ vcf_to_dosage <- function(vcf, ...) { #### Input Validation #### - format <- base::match.arg(format) - marker_names <- base::match.arg(marker_names) + format <- match.arg(format) + marker_names <- match.arg(marker_names) - if (!base::requireNamespace("vcfR", quietly = TRUE)) + if (!requireNamespace("vcfR", quietly = TRUE)) stop("Package 'vcfR' is required to read VCF files. Install it with install.packages(\"vcfR\").") - if (!base::is.null(ploidy)) .check_ploidy(ploidy) - if (!base::is.null(dosage_field) && base::is.null(ploidy)) + if (!is.null(ploidy)) .check_ploidy(ploidy) + if (!is.null(dosage_field) && is.null(ploidy)) stop("ploidy must be supplied when dosage_field is used; it cannot be inferred from a dosage field.") for (arg in c("min_marker_call_rate", "min_sample_call_rate")) { - val <- base::get(arg) - if (base::length(val) != 1 || !base::is.numeric(val) || base::is.na(val) || + val <- get(arg) + if (length(val) != 1 || !is.numeric(val) || is.na(val) || val < 0 || val > 1) stop(arg, " must be a single number between 0 and 1.") } #### Read the VCF #### - if (base::is.character(vcf) && base::length(vcf) == 1) { - if (!base::file.exists(vcf)) + if (is.character(vcf) && length(vcf) == 1) { + if (!file.exists(vcf)) stop("VCF file not found: ", vcf) vcf_obj <- vcfR::read.vcfR(vcf, verbose = verbose, ...) - } else if (base::inherits(vcf, "vcfR")) { + } else if (inherits(vcf, "vcfR")) { vcf_obj <- vcf } else { stop("vcf must be a path to a .vcf / .vcf.gz file or a vcfR object.") @@ -136,23 +136,23 @@ vcf_to_dosage <- function(vcf, # getFIX() drops to a named vector when the VCF has one variant; keep a matrix get_fix <- function(obj) { f <- vcfR::getFIX(obj) - if (!base::is.null(f) && base::is.null(base::dim(f))) - f <- base::matrix(f, nrow = 1, dimnames = base::list(NULL, base::names(f))) + if (!is.null(f) && is.null(dim(f))) + f <- matrix(f, nrow = 1, dimnames = list(NULL, names(f))) f } fix <- get_fix(vcf_obj) - if (base::is.null(fix) || base::nrow(fix) == 0) + if (is.null(fix) || nrow(fix) == 0) stop("The VCF contains no variants.") #### Drop multiallelic markers #### - n_markers_in <- base::nrow(fix) + n_markers_in <- nrow(fix) if (biallelic_only) { - keep <- !base::grepl(",", fix[, "ALT"], fixed = TRUE) & !base::is.na(fix[, "ALT"]) - if (base::sum(keep) == 0) + keep <- !grepl(",", fix[, "ALT"], fixed = TRUE) & !is.na(fix[, "ALT"]) + if (sum(keep) == 0) stop("No biallelic markers remain. Set biallelic_only = FALSE to keep multiallelic markers.") - if (base::any(!keep)) { - warning(base::sum(!keep), " multiallelic marker(s) were dropped. ", + if (any(!keep)) { + warning(sum(!keep), " multiallelic marker(s) were dropped. ", "Set biallelic_only = FALSE to keep them.", call. = FALSE) vcf_obj <- vcf_obj[keep, ] fix <- get_fix(vcf_obj) @@ -161,9 +161,9 @@ vcf_to_dosage <- function(vcf, #### Marker names #### ids <- fix[, "ID"] - positions <- base::paste(fix[, "CHROM"], fix[, "POS"], sep = "_") - id_usable <- !base::any(base::is.na(ids) | ids == "." | ids == "") && - !base::any(base::duplicated(ids)) + positions <- paste(fix[, "CHROM"], fix[, "POS"], sep = "_") + id_usable <- !any(is.na(ids) | ids == "." | ids == "") && + !any(duplicated(ids)) if (marker_names == "id") { if (!id_usable) @@ -174,21 +174,21 @@ vcf_to_dosage <- function(vcf, } else { snp_names <- if (id_usable) ids else positions } - if (base::any(base::duplicated(snp_names))) { + if (any(duplicated(snp_names))) { warning("Duplicated marker names were made unique with a numeric suffix.", call. = FALSE) - snp_names <- base::make.unique(snp_names, sep = "_") + snp_names <- make.unique(snp_names, sep = "_") } #### Extract genotype calls #### - if (!base::is.null(dosage_field)) { + if (!is.null(dosage_field)) { # Dosages are already numeric in the requested FORMAT field dose <- vcfR::extract.gt(vcf_obj, element = dosage_field, as.numeric = TRUE, IDtoRowNames = FALSE) - if (base::all(base::is.na(dose))) + if (all(is.na(dose))) stop("FORMAT field '", dosage_field, "' is absent or empty in this VCF.") - obs_max <- base::suppressWarnings(base::max(dose, na.rm = TRUE)) - if (base::is.finite(obs_max) && obs_max > ploidy) + obs_max <- suppressWarnings(max(dose, na.rm = TRUE)) + if (is.finite(obs_max) && obs_max > ploidy) warning("Dosages above the requested ploidy (", ploidy, ") were found in field '", dosage_field, "'.", call. = FALSE) n_wrong_ploidy <- 0L @@ -196,47 +196,47 @@ vcf_to_dosage <- function(vcf, } else { gt <- vcfR::extract.gt(vcf_obj, element = "GT", IDtoRowNames = FALSE) - if (base::is.null(gt) || base::all(base::is.na(gt))) + if (is.null(gt) || all(is.na(gt))) stop("No GT field found in this VCF.") # Genotype strings repeat heavily, so the parsing is done once per unique # call and mapped back to the full matrix. This keeps large files cheap. - calls <- base::unique(base::as.vector(gt)) - idx <- base::match(gt, calls) - n_obs <- base::tabulate(idx, nbins = base::length(calls)) + calls <- unique(as.vector(gt)) + idx <- match(gt, calls) + n_obs <- tabulate(idx, nbins = length(calls)) - alleles <- base::strsplit(calls, "[/|]") # phased or unphased - n_slots <- base::lengths(alleles) - is_miss <- base::is.na(calls) | n_slots == 0 | - base::vapply(alleles, - function(a) base::any(base::is.na(a) | a == "."), - base::logical(1)) + alleles <- strsplit(calls, "[/|]") # phased or unphased + n_slots <- lengths(alleles) + is_miss <- is.na(calls) | n_slots == 0 | + vapply(alleles, + function(a) any(is.na(a) | a == "."), + logical(1)) # dosage = allele slots that are not the reference allele (0) - dose_u <- base::rep(NA_real_, base::length(calls)) - dose_u[!is_miss] <- base::vapply(alleles[!is_miss], - function(a) base::sum(a != "0"), - base::numeric(1)) + dose_u <- rep(NA_real_, length(calls)) + dose_u[!is_miss] <- vapply(alleles[!is_miss], + function(a) sum(a != "0"), + numeric(1)) #### Ploidy #### obs_slots <- n_slots[!is_miss] obs_count <- n_obs[!is_miss] - if (base::length(obs_slots) == 0) + if (length(obs_slots) == 0) stop("Every genotype call is missing; nothing to convert.") - slot_totals <- base::tapply(obs_count, obs_slots, base::sum) - obs_ploidy <- base::as.integer(base::names(slot_totals)[base::which.max(slot_totals)]) + slot_totals <- tapply(obs_count, obs_slots, sum) + obs_ploidy <- as.integer(names(slot_totals)[which.max(slot_totals)]) ploidy_mismatch <- FALSE - if (base::is.null(ploidy)) { + if (is.null(ploidy)) { ploidy <- obs_ploidy if (verbose) - base::cat("Inferred ploidy from the GT field:", ploidy, "\n") + cat(sprintf("Inferred ploidy from the GT field: %s\n", ploidy)) } else { ploidy_mismatch <- obs_ploidy != ploidy } off_ploidy <- !is_miss & n_slots != ploidy - n_wrong_ploidy <- base::sum(n_obs[off_ploidy]) + n_wrong_ploidy <- sum(n_obs[off_ploidy]) # One warning per situation: a contradicting ploidy, or a few stray calls if (ploidy_mismatch) { @@ -250,58 +250,56 @@ vcf_to_dosage <- function(vcf, } dose_u[is_miss | off_ploidy] <- NA_real_ - dose <- base::matrix(dose_u[idx], - nrow = base::nrow(gt), ncol = base::ncol(gt), - dimnames = base::list(NULL, base::colnames(gt))) - base::storage.mode(dose) <- "integer" # dosages are whole numbers + dose <- matrix(dose_u[idx], + nrow = nrow(gt), ncol = ncol(gt), + dimnames = list(NULL, colnames(gt))) + storage.mode(dose) <- "integer" # dosages are whole numbers } #### Reference instead of alternate dosage #### - ploidy <- base::as.integer(ploidy) + ploidy <- as.integer(ploidy) if (ref_dosage) dose <- ploidy - dose #### Individuals in rows, markers in columns #### - base::rownames(dose) <- snp_names - geno <- base::t(dose) - if (base::is.null(base::rownames(geno)) || base::any(base::rownames(geno) == "")) + rownames(dose) <- snp_names + geno <- t(dose) + if (is.null(rownames(geno)) || any(rownames(geno) == "")) stop("Sample names could not be read from the VCF.") #### Call-rate filtering #### - n_samples_in <- base::nrow(geno) - n_snps_in <- base::ncol(geno) + n_samples_in <- nrow(geno) + n_snps_in <- ncol(geno) if (min_marker_call_rate > 0) { - marker_cr <- base::colMeans(!base::is.na(geno)) + marker_cr <- colMeans(!is.na(geno)) geno <- geno[, marker_cr >= min_marker_call_rate, drop = FALSE] - if (base::ncol(geno) == 0) + if (ncol(geno) == 0) stop("No markers passed min_marker_call_rate = ", min_marker_call_rate, ".") } if (min_sample_call_rate > 0) { - sample_cr <- base::rowMeans(!base::is.na(geno)) + sample_cr <- rowMeans(!is.na(geno)) geno <- geno[sample_cr >= min_sample_call_rate, , drop = FALSE] - if (base::nrow(geno) == 0) + if (nrow(geno) == 0) stop("No individuals passed min_sample_call_rate = ", min_sample_call_rate, ".") } #### Verbose report #### if (verbose) { - base::cat("\n=== VCF Conversion Report ===\n") - base::cat("Ploidy: ", ploidy, "\n") - base::cat("Markers read: ", n_markers_in, "\n") - base::cat("Markers returned: ", base::ncol(geno), "\n") - base::cat("Individuals returned: ", base::nrow(geno), " of ", n_samples_in, "\n", sep = "") - base::cat(base::sprintf("Missing calls: %.2f%%\n", - 100 * base::mean(base::is.na(geno)))) + cat("\n=== VCF Conversion Report ===\n") + cat(sprintf("Ploidy: %s\n", ploidy)) + cat(sprintf("Markers read: %s\n", n_markers_in)) + cat(sprintf("Markers returned: %d\n", ncol(geno))) + cat(sprintf("Individuals returned: %d of %s\n", nrow(geno), n_samples_in)) + cat(sprintf("Missing calls: %.2f%%\n", 100 * mean(is.na(geno)))) if (n_wrong_ploidy > 0) - base::cat("Calls set to missing for wrong ploidy:", n_wrong_ploidy, "\n") - base::cat("Coding: dosage of the ", - if (ref_dosage) "reference" else "alternate", " allele (0..", ploidy, ")\n", - sep = "") + cat(sprintf("Calls set to missing for wrong ploidy: %s\n", n_wrong_ploidy)) + cat(sprintf("Coding: dosage of the %s allele (0..%s)\n", + if (ref_dosage) "reference" else "alternate", ploidy)) } #### Return in the requested layout #### if (format == "matrix") { - base::attr(geno, "ploidy") <- ploidy + attr(geno, "ploidy") <- ploidy return(geno) } diff --git a/man/find_parentage.Rd b/man/find_parentage.Rd index 08e1eae..078bd13 100644 --- a/man/find_parentage.Rd +++ b/man/find_parentage.Rd @@ -86,7 +86,7 @@ have \code{NA} in the remaining ranks. Status refers to rank 1 only.} \item{founders_file}{Optional. Founder IDs as a path to a one-column file, a character vector of ids, or a data.frame (first column). Only used when \code{progeny_file} is a pedigree: founders with -both parents unknown are skipped instead of searched. Default is NULL.} +both parents unknown are skipped and not searched. Default is NULL.} } \value{ A named list (returned invisibly) with elements: diff --git a/man/ped_to_dosage.Rd b/man/ped_to_dosage.Rd index e71c895..61784b9 100644 --- a/man/ped_to_dosage.Rd +++ b/man/ped_to_dosage.Rd @@ -51,7 +51,7 @@ columns, coded as the number of copies of the counted allele (0, 1, 2) with NA for missing calls. PLINK \code{.ped} files are diploid only. } \details{ -The \code{.ped} format stores allele letters rather than a reference/alternate +The \code{.ped} format stores allele letters and has no reference/alternate pair, so the counted allele (allele B) is chosen per marker. By default it is the alphabetically / numerically last allele observed at the marker, so \code{1/2} coded files count allele \code{2} and an \code{A/G} marker counts \code{G}. To code a From f9468fb3e4956c5af28a88a81994c7df8a674815 Mon Sep 17 00:00:00 2001 From: "josue.chinchilla" Date: Mon, 12 Oct 2026 16:25:13 -0400 Subject: [PATCH 8/8] BIGpopA 2.2.0: add user guide in docs, remove stray files --- .DS_Store | Bin 18436 -> 0 bytes R/.DS_Store | Bin 6148 -> 0 bytes _to_delete/test-founders_input.R | 9 - docs/.DS_Store | Bin 6148 -> 0 bytes docs/BIGpopA_user_guide.Rmd | 427 +++++++++++++++++++++++++++++++ docs/BIcolors_template.css | 15 +- inst/.DS_Store | Bin 6148 -> 0 bytes tests/.DS_Store | Bin 8196 -> 0 bytes tests/testthat/.DS_Store | Bin 6148 -> 0 bytes 9 files changed, 439 insertions(+), 12 deletions(-) delete mode 100644 .DS_Store delete mode 100644 R/.DS_Store delete mode 100644 _to_delete/test-founders_input.R delete mode 100644 docs/.DS_Store create mode 100644 docs/BIGpopA_user_guide.Rmd delete mode 100644 inst/.DS_Store delete mode 100644 tests/.DS_Store delete mode 100644 tests/testthat/.DS_Store diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index ee8ada33f742f2eab690d39af43dd298b4197798..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 18436 zcmeHO3v?7$8UFuE2(!!tCXfKTO(2k#X9!84kN|BruR;irkOYzdCE4sQuw}CgyPJeI z#il-LTR@MfwZ8h$qgFkl)`G2lSRd5Z+Ipn43ciaTt4ELRX^)lavHtho*9~CkT2>Dc7?;il(y*#7!~r9`ob9K22W!Du~a-3V&I1;U8B}VkAwA8X&?3s*fU_y zfIS2D4A?Vp{xiU_n+`>lseQF)z@7nn2Fw}Y>q7yZJarrCfK0{GK^Jbm^LlQ=i-YdT z4-jPqM(Q@w0hx-f)JJrsRCYzf#X#Ae$_sdNO5H{}AXC{LP<972-i(Hag5vGeuRt^h zbj#Gf+B0C!K} zUl|NG2X28SK?zlDo1@`KtXe#*EkMtgttRnzH0VWLQq|SkD@Uc*8->erR9Blf+}%w>i5g|sUgvUcY86r>UAa^e{wy4aU%>C+uO#9B z!7-@dR4l-`Sd8U(EnY`5U5&N49vg52Nj8`3ZP7B zCkf~BeGh&LKaHQoFXET*LHr7S3m?X#_-%ZYB>o5Z6MP9z;prqfUY97x_-3LU_cljD z!Ek5?m2$~Kt1oIuIy@4W^q1qnF-}pH8=46l^^h+!6en`MRKP#hilwEq4jWlM_+ zmMmVjuc)BEqO649FWtA#q2v}WTG7<8E7%EH-Tqewwer*lJ-LJ&G!1l=RWG7&)?0x`>= z2G&9+M4%rI!aeW+A?7i78eW3m5L&(gXAm(H$KiNXF$eQ;I;|=*2`iV;y0RQsU=1N9 zUs1LgtBH?RlL$r$D|=}r`2gOIci^3P$XG@0rB&n`_)UBWAHyHxbND>IfIr7q@HGAb z|Ay}{V93U^9G1(buoYgq9s7;oKJN#t`A<9+L4}a3A620NH~aEf3;{@gqXcKZ9J(AH*+{UHCPDmX8Ru z{2?AA`|w$k{1@?8_(%K;zL6x)>5wyxJSQQje29IuXW(6)0Zx}?kqcG_tH=cl(^(0{ z&_rp3ayQlm!(QpY2wtv=S{J_~4)MEn@KQU9+rvV-3sVfgn@qwG57*OEy1|#!|DmM* zDnVH-G($DLu?`xcg{}>-KCu_hO6rvXK6;lIqKWMorv$CP&TRKTVw{@eiJPSSKQIf= zWr(VVA$plnH)^hxx=Dt#^EEsZsheRKo|3vbYj~;)C3SPw@KjR`Pc@z4sb0p0XNr?B z4C>~zLEXG%Qa2M!>Sk^Vb+f^wZoDRS(`!;U`%|c!$IeOJWM2SvV^yXh%}>=ZKmYHE z8}+=*O*PTnREIJ*)kJeswVIo%nqY3K>P4EH=O^ET5jSN9apN|L8%{||;-Ok-JW0-MUFF&Ddn6|+UGoK>aS`(UGdM-@j#uGL;wUqToz!s##XCxH zDaG)v$5kZJ^|%f=9onw366LnL;|=#+2eD;50U;C&DIhXZBu;Gi*Hg z@=EyLKL_omucQg#o_HhORYPs!P=yC>9V*{{Ukk0&K8=@;hms}+p(*is!c(_Ck6Ku4 zjiV!OYfLe0jS_BYEUPFlDyUdeo)kX1L;kUY8X9*}J&hx7VsJU5?755!dqai2p~Bu! zadufwx1mC+k?oI*?MC)Sij=KFgZbm5Wg`X2Z<0abAilenkhVskZW;__+$Xmpq-9f4 zB9;<&4>-YMhFDr&{+Md}|4CacZU29??|r-o|H$_LFWI)?^C#16|3CfC$BU3}+y9q7 zOQ0x ImHz*K0f?L~H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0 `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()` +5. What is the line / breed composition of each sample? -> `allele_freq_poly()`, `solve_composition_poly()`, `assign_composition_line()` + +| Function | Use it to | Main output | +|---|---|---| +| `vcf_to_dosage()` | Convert a VCF to allele-B dosages | Dosage table or matrix | +| `ped_to_dosage()` | Convert a PLINK `.ped` (diploid) to allele-B dosages | Dosage table or matrix | +| `check_ped()` | Fix structural pedigree errors before any genotype work | Corrected pedigree, issue tables | +| `validate_pedigree()` | Test recorded trios against genotypes | Status per trio, corrected pedigree, marker summary | +| `find_parentage()` | Assign or fill in parents from a candidate pool | Best parent(s), error %, ties | +| `allele_freq_poly()` | Allele frequencies of reference lines | Frequency table (markers x lines) | +| `solve_composition_poly()` | Composition of samples from reference frequencies | Proportions and R2 per sample | +| `assign_composition_line()` | Percentages and assigned line per sample | Composition table | + +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. +- **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). + +BIGpopA does not call dosages; use your pipeline first (we recommend [BIGapp](https://github.com/Breeding-Insight/BIGapp)). + +# Install + +BIGpopA requires R >= 4.4.0. Install it from CRAN by copying and pasting this code into your R console: + +```{r, eval=FALSE} +install.packages("BIGpopA") +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. + +| 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()` | + +```{r} +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 +``` + +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`). + +## 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()`. + +```{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) + +# 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: + +| 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) | + +# Pedigree cleaning + +`check_ped()` works on the columns `id`, `male_parent` and `female_parent`. It applies the checks in this order: + +| Issue | Meaning | Action | +|---|---|---| +| Exact duplicates | The identical row appears more than once | Always removed | +| 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 | + +Unknown parents must be coded `0`. + +```{r} +clean_ped_results = check_ped(pedigree, verbose = FALSE) + +# 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. + +```{r} +ped_validate_results = validate_pedigree(clean_ped, genotypes, + verbose = FALSE, plot_results = TRUE) + +ped_validate_report = ped_validate_results$full_results +ped_validate_report +``` + +## Parameters + +| Argument | Default | Effect | +|---|---|---| +| `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 | + +## Status of each trio + +| Status | Meaning | +|---|---| +| `pass` | Mendelian error at or below `trio_error_threshold` | +| `fail` | 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 | + +In the toy data `C1`, `C2` and `C3` pass, `C_bad` fails, `C_low` is `low_markers` and `C_missing` has `no_genotype_data`. + +## 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: + +| `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 | +| `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) | + +```{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. + +## Marker summary + +`marker_summary` lists every marker with `trios_tested`, `trios_mismatch` and `mismatch_pct`, sorted worst first. It answers a different question than the trio table: *is the problem the individual or the 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. + +```{r} +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`. + +```{r} +mt_results = validate_pedigree(clean_ped, genotypes, verbose = FALSE, + plot_results = FALSE, marker_trio_table = TRUE) +mt_results$marker_trio_table +``` + +| Label | Meaning | +|---|---| +| `match` | The trio is consistent at this marker | +| `mismatch_male`, `mismatch_female` | Only that parent is inconsistent with the offspring | +| `mismatch_both` | Each parent alone is inconsistent with the offspring | +| `mismatch_combination` | Each parent alone is compatible, but the pair cannot produce the offspring | +| `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 + +Missing marker calls (`NA`) are handled marker by marker and never count as errors: + +- **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. + +# Parentage assignment + +When parentage is unknown, `find_parentage()` searches a pool of candidate parents for the best assignment to each progeny. + +```{r} +find_parentage_results = find_parentage(genotypes, parents, progeny, + method = "best_pair", + verbose = FALSE, plot_results = TRUE) + +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`. + +## 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`). + +```{r} +find_parentage(genotypes, parents, progeny, method = "best_pair", top_n = 3, + verbose = FALSE, plot_results = FALSE)$full_results +``` + +## 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`. + +```{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 +``` + +# 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. + +```{r} +# 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 +``` + +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. + +```{r} +# 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. + +`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. + +```{r} +pred_results = assign_composition_line(solve_composition_poly(validation, freq, ploidy = 4)) +pred_results +``` + +`Sample1` is assigned Jewel and `Sample2` Draper (each essentially 100%), while the admixed `Sample3` splits roughly evenly between them. For diploid data set `ploidy = 2`. Save the table with `write.table(pred_results, "composition.txt", sep = "\t", row.names = FALSE)`. + +# Complete workflow + +```{r, eval=FALSE} +# genotypes +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) + +# breed/line composition +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) +``` + +# Tips + +- 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. +- 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. + +# How to cite + +## BIGpopA + +Chinchilla-Vargas J, Sandercock AM (2026). *BIGpopA: Pedigree Validation Genetic Composition of Diploids & Polyploids*. R package version 2.2.0. CRAN: + +Run `citation("BIGpopA")` in R for the BibTeX entry. + +## Methods + +- Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP (2017). Estimation of genome-wide and locus-specific breed composition in pigs. *Translational Animal Science* 1(1):36-44. + +# Session information + +```{r} +sessionInfo() +``` diff --git a/docs/BIcolors_template.css b/docs/BIcolors_template.css index 0d12431..a6fdb41 100644 --- a/docs/BIcolors_template.css +++ b/docs/BIcolors_template.css @@ -280,8 +280,8 @@ pre code { /* Inline code styling */ p code, li code, td code, span code { - background-color: var(--code-bg-dark) !important; - color: white !important; + background-color: #eef2f2 !important; + color: #1a1a1a !important; padding: 2px 6px !important; border-radius: 4px !important; font-family: 'Fira Code', 'Monaco', 'Consolas', monospace !important; @@ -485,7 +485,7 @@ body::after { url('usda-logo-color.png'), url('IFAS.jpg'), url('cornell_seal_simple_web_b31b1b.png'); - background-size: 150px 80px, 150px 80px, 100px 80px; + background-size: auto 55px, 150px 80px, 100px 80px; background-repeat: no-repeat, no-repeat, no-repeat; background-position: calc(100% - 530px) center, @@ -493,3 +493,12 @@ body::after { calc(100% - 240px) center; box-sizing: border-box; } + +/* Table borders for readability */ +.table, table.table { + border: 1px solid #b0b0b0; + border-collapse: collapse; +} +.table th, .table td, table.table th, table.table td { + border: 1px solid #b0b0b0 !important; +} diff --git a/inst/.DS_Store b/inst/.DS_Store deleted file mode 100644 index 5db1497c40e49c85f0efe88eb2f3feec7d8c6fde..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeHKO-sW-5Pe&FsCenoV}5}_`~$H@YC$Sk^{Q4w{iwFI>A_omx^H$yX-P%B6pNPl>$zMbmZ47gtlN6hj9htwX<&|k_b-%wUv zW?tfiObnT$@_o*k@FAQ^-s{f!k!S50{au^^XTTZwJp-)SBCQ=oADsbbz!~^tK)w$N zRWLKG6yv9ZLu>(v@*5e2bLk}{CNazmD@8t`1Pvu>Xs{L|XgKYO`DKQcqJ|?_^AYUL z!8()>-_HG$=#G#n`sfTe1FH;d>9Q^L|M=_ve>KQg&VV!UuNVl;e2{mTlGWB)a#Cv( r>J3#y{3^xI5Kdw#rmvLZU1}8C6PXY*!%C4BivI{i8hmgDHp;*|5e{Ju diff --git a/tests/.DS_Store b/tests/.DS_Store deleted file mode 100644 index 7055d5c12008d702ecacc98e9b571cf77986776c..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 8196 zcmeHMTWl0n7(U;$g&8N%X(<=SEv{6lKm$d&v?#FMwjva0?Y8v7rOfUObi#C|?#yn% zidLUg)JHWMNC_H`aRlNB%!~l9-DRrVE#}jBS^b`Mye!GiWJev_4;mZ4fmBvrF=wtkS61cy znPG1%6Bt3kPaFBYeAVN5PPQ~RvU_Z6Y)GqV_1wU)T`M0N*rr7)`}%Cx^v2R&!8QHR zHX{VGA}d4MqVe%Nnm05g8;>_6$2T@NH6)wYH61^$$P4S%+|@mF#2InD$Alw@J_c;w zG{+vF%`4N-Iazv2ZZyMX$1k=Dm&?|dT)t3pS-C+~=jlDYd;0nZhBU3#dA=Rko;%PV zSVsdc;~Py8^$+@9K4=X$>$hotp;gV#=p5z7{vIpu6IGOx^A{&>xozc|bAndHbW&*ZzA+G$znFI%eXd%2`FZ_#K*2FvSpwZ|v7G$T@6 z)u5_9N9gTqesomhT2=3LjI8BorXtQ+zd_aeioz?cBVul|-YsiqWqRwjTeJ!B6l(!|xvm7dy3|TBt@|jO3TI(J>r#z^T_eck$nUyZ*2vOq z2Rpz<*)ev4J;PpNr`TzBhJC`mWaruU>?d}S{m%Yim)M^u$2?R*Lltz?V+B@X71kkz z7NpUEE_5S<9_&XJCM+C+hXTfM6vyx=p1_lM3Qyx%yo^`yD&E3-cpo3&Lwt&}_zLIn z9e%;DxQI(qxl|#|kt!uks+JZ=%cXj$QEHab(l)6}JiS!PxWP$p!ii)MnJuYCMN&y8 zXHTlcwtKhl=-7G9Ce``*B8?Ue|Ge z3MHUcCYB9}jf$*6EH7P5!5F{y96#zGXkKp9$xG zp$sw>BZ=Fw8co=Ot%U4Wv|~3qaUbEEqx%33!axoVMhM}f$YUH2;vqbYM+oI72m|=ylX0&@!M=@E@tws@BjGzFIC0=t^5CPYgKSc diff --git a/tests/testthat/.DS_Store b/tests/testthat/.DS_Store deleted file mode 100644 index 5008ddfcf53c02e82d7eee2e57c38e5672ef89f6..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 6148 zcmeH~Jr2S!425mzP>H1@V-^m;4Wg<&0T*E43hX&L&p$$qDprKhvt+--jT7}7np#A3 zem<@ulZcFPQ@L2!n>{z**++&mCkOWA81W14cNZlEfg7;MkzE(HCqgga^y>{tEnwC%0;vJ&^%eQ zLs35+`xjp>T0