diff --git a/.DS_Store b/.DS_Store deleted file mode 100644 index ee8ada3..0000000 Binary files a/.DS_Store and /dev/null differ diff --git a/DESCRIPTION b/DESCRIPTION index 5e9c25b..c55da55 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.2.0 Authors@R: c(person(given='Josue', family='Chinchilla-Vargas', email='josue.chinchilla@ufl.edu', 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 f8366fc..56ae501 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,18 @@ +# 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"`. +- `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. +- `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/.DS_Store b/R/.DS_Store deleted file mode 100644 index 5008ddf..0000000 Binary files a/R/.DS_Store and /dev/null differ 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/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/find_parentage.R b/R/find_parentage.R index 4d1ab0f..fc09634 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 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 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 +#' `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 != round(top_n)) + stop("top_n must be a single positive integer.") + top_n <- as.integer(top_n) # Accept file path OR in-memory data.frame / data.table read_flex <- function(x, label, ...) { @@ -143,107 +185,170 @@ 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) is.na(x) | trimws(as.character(x)) %in% c("", "0", "NA") + method_reported <- method + pedigree_mode <- method == "fill_pedigree" + if (pedigree_mode) { + 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 := 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) + 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 (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_ + status = NA_character_, + tied_candidates = NA_integer_ ) + 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_ } - best_idx <- base::which.min(percent_mismatch) - if (base::length(best_idx) == 0) { + if (top_n > 1) { + valid_idx <- which(!is.na(percent_mismatch)) + ranked_idx <- valid_idx[order(percent_mismatch[valid_idx], + -comparisons[valid_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 <- which.min(percent_mismatch) + } + 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", + as.integer(sum(percent_mismatch == percent_mismatch[best_idx], + na.rm = TRUE))) + + 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(Filter(Negate(is.null), extra_list), + fill = TRUE, use.names = TRUE) + if (top_n > 1) { + for (k in seq(2, top_n)) { + cols <- paste0(c("best_match_", "mendelian_error_pct_", "markers_tested_"), k) + for (col in cols) { + 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) + } + } } final_df <- results_dt } @@ -254,25 +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 - progeny_mat <- genos_mat[progeny_ids, , drop = FALSE] - n_progeny <- base::nrow(progeny_mat) - n_pairs <- base::nrow(parent_pairs) + n_all <- length(progeny_ids) + if (pedigree_mode) { + kmale <- ifelse(is_unknown(progeny_candidates$male_parent), + NA_character_, progeny_candidates$male_parent) + kfem <- ifelse(is_unknown(progeny_candidates$female_parent), + NA_character_, progeny_candidates$female_parent) + } else { + kmale <- kfem <- rep(NA_character_, n_all) + } + 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 <- 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 @@ -281,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) ) @@ -297,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, @@ -305,76 +423,151 @@ 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 <- ifelse(search_mode == "known_male", kmale, + ifelse(search_mode == "known_female", kfem, + NA_character_)) + if (pedigree_mode) { + known_err <- rep(NA_real_, n_all) + for (j in which(search_mode != "pair")) { + kid <- known_id_vec[j] + if (kid %in% rownames(genos_hom_mat)) { + kh <- genos_hom_mat[kid, ] + ph <- genos_hom_mat[progeny_ids[j], ] + n_cmp <- sum(!is.na(kh) & !is.na(ph)) + known_err[j] <- if (n_cmp == 0) NA_real_ else + round(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 <- 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 <- 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% 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% 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 + } + 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 <- 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[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_ } - 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 } - 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 <- which(!is.na(percent_mismatch)) + ranked_idx <- valid_idx[order(percent_mismatch[valid_idx], + -comparisons[valid_idx])] + best_indices <- ranked_idx[seq_len(min(top_n, length(ranked_idx)))] + } else { + best_indices <- which(percent_mismatch == min_mismatch_val) + if (length(best_indices) > 1) { + best_markers_per_pair <- comparisons[best_indices] + max_markers <- 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) + best_error <- round(min_mismatch_val, 2) a_status <- assign_status(best_markers, best_error) - if (!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) } - num_to_report <- base::min(base::nrow(best_pairs), - if (show_ties) base::nrow(best_pairs) else 1) + report_all <- show_ties || top_n > 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", + as.integer(sum(percent_mismatch == min_mismatch_val, + na.rm = TRUE))) - if (show_ties && 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("markers_tested_", k)]] <- comparisons[best_indices[k]] + if (report_all && num_to_report > 1) { + 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 @@ -392,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, " | ", - "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 #### @@ -435,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, " | ", @@ -476,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", @@ -485,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 548241a..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 } @@ -301,3 +301,73 @@ 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 >= pmax(0, g - h) & offspring <= 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 <- 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", + !mm ~ "match", + male_bad & !female_bad ~ "mismatch_male", + female_bad & !male_bad ~ "mismatch_female", + male_bad & female_bad ~ "mismatch_both", + 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]])) + } + # 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 d0e7046..56805ad 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. @@ -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. @@ -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{ @@ -47,9 +50,48 @@ #' \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{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 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.} +#' \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( @@ -95,11 +137,13 @@ 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 plot_status <- recommended_correction <- NULL + marker <- trios_tested <- trios_mismatch <- mismatch_pct <- NULL #### Input validation #### if (trio_error_threshold < 0 || trio_error_threshold > 100) @@ -134,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)] @@ -146,24 +190,15 @@ 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") + 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) @@ -173,35 +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.") - #### 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 <- 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] @@ -211,10 +231,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) { @@ -224,37 +240,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 <- 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 + + as.numeric((mm_vec %in% TRUE) & tt_vec) if (markers_tested == 0) { status <- "no_data" correction_decision <- "none" @@ -272,40 +276,32 @@ 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 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, 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" - 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) @@ -318,23 +314,19 @@ 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, 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 ) }) 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, @@ -344,60 +336,44 @@ 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)) } #### 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) - if (decision == "remove_male_parent") { - data.table::set(corrected_pedigree, row_idx, "male_parent", "0") - } else if (decision == "remove_female_parent") { + row_idx <- which(corrected_pedigree$id == prog_id) + # low_markers_* decisions apply the same removal as the plain 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")) 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 #### 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 #### @@ -406,18 +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)] - plot_df$plot_status <- dplyr::case_when( - plot_df$recommended_correction %in% c("none", "keep_both", - "low_markers_keep_both") ~ "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", - TRUE ~ "other" - ) + # Colour by the trio error first, so no trio above the threshold is shown as pass + 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) @@ -438,12 +417,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, @@ -458,8 +440,34 @@ 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 <- 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)) + 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( + c(list(marker = marker_names), status_list)) + } + + #### Per-marker mismatch summary #### + marker_summary <- data.table::data.table( + marker = marker_names, + 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"], @@ -470,7 +478,9 @@ validate_pedigree <- function(pedigree_file, genotypes_file, "missing_female_parent")], 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)) + 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/README.md b/README.md index ace0952..d84f990 100644 --- a/README.md +++ b/README.md @@ -41,10 +41,11 @@ 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 | +| `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 256213c..ca397e9 100644 --- a/cran-comments.md +++ b/cran-comments.md @@ -2,24 +2,37 @@ 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 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. @@ -32,5 +45,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/docs/.DS_Store b/docs/.DS_Store deleted file mode 100644 index 5008ddf..0000000 Binary files a/docs/.DS_Store and /dev/null differ diff --git a/docs/BIGpopA_user_guide.Rmd b/docs/BIGpopA_user_guide.Rmd new file mode 100644 index 0000000..291a5ec --- /dev/null +++ b/docs/BIGpopA_user_guide.Rmd @@ -0,0 +1,427 @@ +--- +title: "BIGpopA Tutorial" +subtitle: "Pedigree cleaning, validation, parentage assignment and breed composition from SNP genotype data" +author: "Josue Chinchilla-Vargas and Alexander Sandercock" +date: "`r Sys.Date()`" +output: + html_document: + theme: flatly + highlight: tango + css: BIcolors_template.css + includes: + in_header: header.html + toc: true + toc_depth: 3 + toc_float: + collapsed: true + smooth_scroll: true + number_sections: true + code_folding: show + df_print: paged +--- + +```{r, echo=FALSE} +knitr::opts_chunk$set( + echo = TRUE, message = FALSE, warning = FALSE, + fig.align = "center", fig.width = 7, fig.height = 4.5 +) +``` + +# Introduction + +BIGpopA is an R package for pedigree quality control, parentage assignment and breed/line composition from SNP genotypes. It works for diploid and polyploid species and is the analytical engine behind the [Familia](https://github.com/Breeding-Insight/Familia) Shiny app. This guide runs the functions in the order a breeding program usually needs them: + +1. Convert genotypes to the dosage table used by every function. -> `vcf_to_dosage()`, `ped_to_dosage()` +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 5db1497..0000000 Binary files a/inst/.DS_Store and /dev/null differ 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 c61ffef..078bd13 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,22 @@ 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{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.} + +\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 and not searched. Default is NULL.} } \value{ A named list (returned invisibly) with elements: @@ -73,7 +94,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.} } } @@ -82,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/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 diff --git a/man/validate_pedigree.Rd b/man/validate_pedigree.Rd index 0954eae..1d0b2bc 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{ @@ -27,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 @@ -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: @@ -66,17 +71,57 @@ 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{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.} } } \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. } +\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{ +\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( diff --git a/tests/.DS_Store b/tests/.DS_Store deleted file mode 100644 index 7055d5c..0000000 Binary files a/tests/.DS_Store and /dev/null differ diff --git a/tests/testthat/.DS_Store b/tests/testthat/.DS_Store deleted file mode 100644 index 5008ddf..0000000 Binary files a/tests/testthat/.DS_Store and /dev/null differ 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-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-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") +}) 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) +}) diff --git a/tests/testthat/test-validate_pedigree.R b/tests/testthat/test-validate_pedigree.R index 7930ae2..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) }) @@ -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))) }) @@ -262,7 +260,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", { @@ -346,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", { @@ -359,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", { @@ -372,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", { @@ -423,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", { @@ -504,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", { @@ -617,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) }) @@ -747,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_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/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 4f81ffc..da847ff 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 @@ -238,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 ``` @@ -272,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