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/NEWS.md b/NEWS.md index 66f583c..56ae501 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,4 +1,5 @@ # BIGpopA 2.2.0 +- `validate_pedigree()` and `find_parentage()`: `founders_file` now also accepts a character vector of ids or a data.frame (first column), in addition to a file path. All inputs can now be supplied from memory. - `validate_pedigree()`: removed the suggested replacement parents (`best_male_candidate`, `best_male_candidate_error_pct`, `best_female_candidate`, `best_female_candidate_error_pct`) from `full_results` and the per-status tables. Use `find_parentage(method = "fill_pedigree")` on `corrected_pedigree` to fill in missing parents. - `validate_pedigree()`: new `marker_summary` element in the returned list gives, for every marker, the number of trios tested, the number with a Mendelian mismatch and the mismatch percentage, sorted worst first. It covers all trios with both parents known and genotyped, and helps spot markers (e.g. paralogous or miscalled) that cause many failures. - `find_parentage()`: new `method = "fill_pedigree"` fills in missing parents of a pedigree. `progeny_file` must have `id`, `male_parent` and `female_parent` columns (e.g. `validate_pedigree()$corrected_pedigree`); parents coded `0`/`NA`/empty are unknown. A known male searches only females, a known female only males, and rows with both unknown get the full pair search (skip founders with the new `founders_file`). Rows with both parents known are skipped. Results gain `search_mode` and `known_parent_error_pct`, and a known parent without genotypes gets `status = "known_parent_no_genotype"`. diff --git a/R/.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/find_parentage.R b/R/find_parentage.R index e679ff9..fc09634 100644 --- a/R/find_parentage.R +++ b/R/find_parentage.R @@ -42,9 +42,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 founders_file Character, optional. Path to a one-column file of -#' founder IDs. Only used when `progeny_file` is a pedigree: founders with -#' both parents unknown are skipped instead of searched. Default is NULL. +#' @param 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 @@ -156,9 +156,9 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, stop("error_threshold must be between 0 and 100.") .check_ploidy(ploidy) if (!is.numeric(top_n) || length(top_n) != 1 || is.na(top_n) || - top_n < 1 || top_n != base::round(top_n)) + top_n < 1 || top_n != round(top_n)) stop("top_n must be a single positive integer.") - top_n <- base::as.integer(top_n) + top_n <- as.integer(top_n) # Accept file path OR in-memory data.frame / data.table read_flex <- function(x, label, ...) { @@ -189,167 +189,163 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, c("id", "male_parent", "female_parent")) #### Pedigree ("fill in the blanks") mode #### - is_unknown <- function(x) base::is.na(x) | base::trimws(base::as.character(x)) %in% c("", "0", "NA") + is_unknown <- function(x) is.na(x) | trimws(as.character(x)) %in% c("", "0", "NA") method_reported <- method pedigree_mode <- method == "fill_pedigree" if (pedigree_mode) { - if (!base::all(c("male_parent", "female_parent") %in% base::names(progeny_candidates))) + if (!all(c("male_parent", "female_parent") %in% names(progeny_candidates))) stop("method = 'fill_pedigree' requires a progeny_file with id, male_parent ", "and female_parent columns.") method <- "best_pair" - progeny_candidates[, male_parent := base::as.character(male_parent)] - progeny_candidates[, female_parent := base::as.character(female_parent)] - founder_ids <- if (!is.null(founders_file)) - base::unique(data.table::fread(founders_file, header = FALSE, - colClasses = "character")[[1]]) else - base::character(0) + progeny_candidates[, male_parent := as.character(male_parent)] + progeny_candidates[, female_parent := as.character(female_parent)] + founder_ids <- .read_founders(founders_file) um <- is_unknown(progeny_candidates$male_parent) uf <- is_unknown(progeny_candidates$female_parent) is_founder <- um & uf & progeny_candidates$id %in% founder_ids if (verbose) - base::cat("Pedigree mode:", base::sum(!um & !uf), "complete trios skipped,", - base::sum(is_founder), "founders skipped,", - base::sum((um | uf) & !is_founder), "rows with a missing parent to search.\n") + cat(sprintf("Pedigree mode: %d complete trios skipped, %d founders skipped, %d rows with a missing parent to search.\n", + sum(!um & !uf), sum(is_founder), sum((um | uf) & !is_founder))) progeny_candidates <- progeny_candidates[(um | uf) & !is_founder] - if (base::nrow(progeny_candidates) == 0) + if (nrow(progeny_candidates) == 0) stop("No progeny with a missing parent to fill in.") } valid_ids <- genos$id - removed_parents <- base::setdiff(all_parents$id, valid_ids) - if (base::length(removed_parents) > 0) { + removed_parents <- setdiff(all_parents$id, valid_ids) + if (length(removed_parents) > 0) { warning("The following parent IDs were not in the genotype file and will not be analyzed: ", paste(removed_parents, collapse = ", "), call. = FALSE) all_parents <- all_parents[id %in% valid_ids] } - removed_progeny <- base::setdiff(progeny_candidates$id, valid_ids) - if (base::length(removed_progeny) > 0) { + removed_progeny <- setdiff(progeny_candidates$id, valid_ids) + if (length(removed_progeny) > 0) { warning("The following progeny IDs were not in the genotype file and will not be analyzed: ", paste(removed_progeny, collapse = ", "), call. = FALSE) progeny_candidates <- progeny_candidates[id %in% valid_ids] } - if (!"sex" %in% base::colnames(all_parents)) { + if (!"sex" %in% colnames(all_parents)) { warning("No 'sex' column in parents file. All parents treated as ambiguous ('A').") all_parents[, sex := "A"] } - all_parents[, sex := base::toupper(sex)] - all_parents <- base::unique(all_parents, by = c("id", "sex")) - male_parent_candidates <- base::unique(all_parents[sex %in% c("M", "A", "NA"), .SD], + all_parents[, sex := toupper(sex)] + all_parents <- unique(all_parents, by = c("id", "sex")) + male_parent_candidates <- unique(all_parents[sex %in% c("M", "A", "NA"), .SD], by = "id") - female_parent_candidates <- base::unique(all_parents[sex %in% c("F", "A", "NA")], + female_parent_candidates <- unique(all_parents[sex %in% c("F", "A", "NA")], by = "id") - if (base::nrow(male_parent_candidates) == 0 && + if (nrow(male_parent_candidates) == 0 && method %in% c("best_male_parent", "best_pair")) warning("No valid male parent candidates remain after filtering.", call. = FALSE) - if (base::nrow(female_parent_candidates) == 0 && + if (nrow(female_parent_candidates) == 0 && method %in% c("best_female_parent", "best_pair")) warning("No valid female parent candidates remain after filtering.", call. = FALSE) - if (base::nrow(progeny_candidates) == 0) + if (nrow(progeny_candidates) == 0) stop("No valid progeny candidates remain after filtering.") #### Pre-compute genotype matrices once #### - genos_mat <- base::as.matrix(genos, rownames = "id") + genos_mat <- as.matrix(genos, rownames = "id") genos_hom <- data.table::copy(genos) - marker_cols <- base::setdiff(base::names(genos_hom), "id") + marker_cols <- setdiff(names(genos_hom), "id") for (col in marker_cols) - genos_hom[base::get(col) > 0 & base::get(col) < ploidy, (col) := NA_integer_] - genos_hom_mat <- base::as.matrix(genos_hom, rownames = "id") + genos_hom[get(col) > 0 & get(col) < ploidy, (col) := NA_integer_] + genos_hom_mat <- as.matrix(genos_hom, rownames = "id") #### Status helper #### assign_status <- function(markers, error_pct) { - base::ifelse(markers < min_markers, "low_markers", - base::ifelse(error_pct > error_threshold, "high_error", "pass")) + ifelse(markers < min_markers, "low_markers", + ifelse(error_pct > error_threshold, "high_error", "pass")) } #### Logic for Homozygous Matching Methods #### if (method %in% c("best_male_parent", "best_female_parent", "best_match")) { - parent_ids <- base::switch(method, + parent_ids <- switch(method, "best_male_parent" = male_parent_candidates$id, "best_female_parent" = female_parent_candidates$id, - "best_match" = base::union(male_parent_candidates$id, + "best_match" = union(male_parent_candidates$id, female_parent_candidates$id)) - parent_genos <- genos_hom_mat[base::rownames(genos_hom_mat) %in% parent_ids, , drop = FALSE] - progeny_genos <- genos_hom_mat[base::rownames(genos_hom_mat) %in% progeny_candidates$id, , drop = FALSE] + parent_genos <- genos_hom_mat[rownames(genos_hom_mat) %in% parent_ids, , drop = FALSE] + progeny_genos <- genos_hom_mat[rownames(genos_hom_mat) %in% progeny_candidates$id, , drop = FALSE] - n_progeny <- base::nrow(progeny_genos) + n_progeny <- nrow(progeny_genos) results_dt <- data.table::data.table( - id = base::rownames(progeny_genos), + id = rownames(progeny_genos), best_match = NA_character_, mendelian_error_pct = NA_real_, markers_tested = NA_integer_, status = NA_character_, tied_candidates = NA_integer_ ) - extra_list <- base::vector("list", n_progeny) + extra_list <- vector("list", n_progeny) - for (i in base::seq_len(n_progeny)) { - prog_id <- base::rownames(progeny_genos)[i] + for (i in seq_len(n_progeny)) { + prog_id <- rownames(progeny_genos)[i] progeny_vec <- progeny_genos[i, ] - progeny_mat <- base::matrix(progeny_vec, - nrow = base::nrow(parent_genos), - ncol = base::ncol(parent_genos), + progeny_mat <- matrix(progeny_vec, + nrow = nrow(parent_genos), + ncol = ncol(parent_genos), byrow = TRUE) - mismatches <- base::rowSums(parent_genos != progeny_mat, na.rm = TRUE) - comparisons <- base::rowSums(!base::is.na(parent_genos) & !base::is.na(progeny_mat)) + mismatches <- rowSums(parent_genos != progeny_mat, na.rm = TRUE) + comparisons <- rowSums(!is.na(parent_genos) & !is.na(progeny_mat)) percent_mismatch <- (mismatches / comparisons) * 100 - percent_mismatch[base::is.nan(percent_mismatch)] <- NA + percent_mismatch[is.nan(percent_mismatch)] <- NA if (exclude_self_match) { - self_idx <- base::rownames(parent_genos) == prog_id + self_idx <- rownames(parent_genos) == prog_id percent_mismatch[self_idx] <- NA_real_ } if (top_n > 1) { - valid_idx <- base::which(!base::is.na(percent_mismatch)) - ranked_idx <- valid_idx[base::order(percent_mismatch[valid_idx], + valid_idx <- which(!is.na(percent_mismatch)) + ranked_idx <- valid_idx[order(percent_mismatch[valid_idx], -comparisons[valid_idx])] - top_idx <- ranked_idx[base::seq_len(base::min(top_n, base::length(ranked_idx)))] - best_idx <- top_idx[base::seq_len(base::min(1L, base::length(top_idx)))] + top_idx <- ranked_idx[seq_len(min(top_n, length(ranked_idx)))] + best_idx <- top_idx[seq_len(min(1L, length(top_idx)))] } else { - best_idx <- base::which.min(percent_mismatch) + best_idx <- which.min(percent_mismatch) } - if (base::length(best_idx) == 0) { + if (length(best_idx) == 0) { data.table::set(results_dt, i, "markers_tested", 0L) data.table::set(results_dt, i, "status", "low_markers") next } best_markers <- comparisons[best_idx] - best_error <- base::round(percent_mismatch[best_idx], 2) - data.table::set(results_dt, i, "best_match", base::rownames(parent_genos)[best_idx]) + best_error <- round(percent_mismatch[best_idx], 2) + data.table::set(results_dt, i, "best_match", rownames(parent_genos)[best_idx]) data.table::set(results_dt, i, "mendelian_error_pct", best_error) - data.table::set(results_dt, i, "markers_tested", base::as.integer(best_markers)) + data.table::set(results_dt, i, "markers_tested", as.integer(best_markers)) data.table::set(results_dt, i, "status", assign_status(best_markers, best_error)) data.table::set(results_dt, i, "tied_candidates", - base::as.integer(base::sum(percent_mismatch == percent_mismatch[best_idx], + as.integer(sum(percent_mismatch == percent_mismatch[best_idx], na.rm = TRUE))) - if (top_n > 1 && base::length(top_idx) > 1) { - rank_row <- base::list(id = prog_id) - for (k in base::seq(2, base::length(top_idx))) { - rank_row[[base::paste0("best_match_", k)]] <- base::rownames(parent_genos)[top_idx[k]] - rank_row[[base::paste0("mendelian_error_pct_", k)]] <- base::round(percent_mismatch[top_idx[k]], 2) - rank_row[[base::paste0("markers_tested_", k)]] <- base::as.integer(comparisons[top_idx[k]]) + if (top_n > 1 && length(top_idx) > 1) { + rank_row <- list(id = prog_id) + for (k in seq(2, length(top_idx))) { + rank_row[[paste0("best_match_", k)]] <- rownames(parent_genos)[top_idx[k]] + rank_row[[paste0("mendelian_error_pct_", k)]] <- round(percent_mismatch[top_idx[k]], 2) + rank_row[[paste0("markers_tested_", k)]] <- as.integer(comparisons[top_idx[k]]) } extra_list[[i]] <- data.table::as.data.table(rank_row) } } # Add ranks 2..top_n as suffix columns, keeping the progeny order - extra_dt <- data.table::rbindlist(base::Filter(Negate(base::is.null), extra_list), + extra_dt <- data.table::rbindlist(Filter(Negate(is.null), extra_list), fill = TRUE, use.names = TRUE) if (top_n > 1) { - for (k in base::seq(2, top_n)) { - cols <- base::paste0(c("best_match_", "mendelian_error_pct_", "markers_tested_"), k) + for (k in seq(2, top_n)) { + cols <- paste0(c("best_match_", "mendelian_error_pct_", "markers_tested_"), k) for (col in cols) { - vals <- if (base::nrow(extra_dt) > 0 && col %in% base::names(extra_dt)) - extra_dt[[col]][base::match(results_dt$id, extra_dt$id)] else - base::rep(if (grepl("^best_match", col)) NA_character_ else - if (grepl("^markers", col)) NA_integer_ else NA_real_, base::nrow(results_dt)) + vals <- if (nrow(extra_dt) > 0 && col %in% names(extra_dt)) + extra_dt[[col]][match(results_dt$id, extra_dt$id)] else + rep(if (grepl("^best_match", col)) NA_character_ else + if (grepl("^markers", col)) NA_integer_ else NA_real_, nrow(results_dt)) data.table::set(results_dt, j = col, value = vals) } } @@ -363,38 +359,38 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, female_parent = female_parent_candidates$id) if (!allow_parent_selfing) { parent_pairs <- parent_pairs[male_parent != female_parent] - if (verbose) base::cat("Parent selfing is disallowed. Pairs with identical parents are removed.\n") + if (verbose) cat("Parent selfing is disallowed. Pairs with identical parents are removed.\n") } - if (base::nrow(parent_pairs) == 0) stop("No valid parent pairs to test.") + if (nrow(parent_pairs) == 0) stop("No valid parent pairs to test.") male_parent_genos_mat <- genos_mat[parent_pairs$male_parent, , drop = FALSE] female_parent_genos_mat <- genos_mat[parent_pairs$female_parent, , drop = FALSE] progeny_ids <- progeny_candidates$id - n_all <- base::length(progeny_ids) + n_all <- length(progeny_ids) if (pedigree_mode) { - kmale <- base::ifelse(is_unknown(progeny_candidates$male_parent), + kmale <- ifelse(is_unknown(progeny_candidates$male_parent), NA_character_, progeny_candidates$male_parent) - kfem <- base::ifelse(is_unknown(progeny_candidates$female_parent), + kfem <- ifelse(is_unknown(progeny_candidates$female_parent), NA_character_, progeny_candidates$female_parent) } else { - kmale <- kfem <- base::rep(NA_character_, n_all) + kmale <- kfem <- rep(NA_character_, n_all) } - search_mode <- base::ifelse(!base::is.na(kmale), "known_male", - base::ifelse(!base::is.na(kfem), "known_female", "pair")) - pair_idx <- base::which(search_mode == "pair") + search_mode <- ifelse(!is.na(kmale), "known_male", + ifelse(!is.na(kfem), "known_female", "pair")) + pair_idx <- which(search_mode == "pair") # Full pair matrices are only built for progeny with no known parent progeny_mat <- genos_mat[progeny_ids[pair_idx], , drop = FALSE] - n_progeny <- base::nrow(progeny_mat) - n_pairs <- base::nrow(parent_pairs) + n_progeny <- nrow(progeny_mat) + n_pairs <- nrow(parent_pairs) - mismatch_mat <- base::matrix( - base::vapply(base::seq_len(n_progeny), function(j) { + mismatch_mat <- matrix( + vapply(seq_len(n_progeny), function(j) { progeny_vec <- progeny_mat[j, ] - progeny_pair_mat <- base::matrix(progeny_vec, + progeny_pair_mat <- matrix(progeny_vec, nrow = n_pairs, - ncol = base::ncol(progeny_mat), + ncol = ncol(progeny_mat), byrow = TRUE) - base::rowSums( + rowSums( .mend_mismatch(male_parent_genos_mat, female_parent_genos_mat, progeny_pair_mat, ploidy), na.rm = TRUE @@ -403,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) ) @@ -419,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, @@ -432,19 +428,19 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, ) # Known-parent bookkeeping (pedigree mode only) - known_id_vec <- base::ifelse(search_mode == "known_male", kmale, - base::ifelse(search_mode == "known_female", kfem, + known_id_vec <- ifelse(search_mode == "known_male", kmale, + ifelse(search_mode == "known_female", kfem, NA_character_)) if (pedigree_mode) { - known_err <- base::rep(NA_real_, n_all) - for (j in base::which(search_mode != "pair")) { + known_err <- rep(NA_real_, n_all) + for (j in which(search_mode != "pair")) { kid <- known_id_vec[j] - if (kid %in% base::rownames(genos_hom_mat)) { + if (kid %in% rownames(genos_hom_mat)) { kh <- genos_hom_mat[kid, ] ph <- genos_hom_mat[progeny_ids[j], ] - n_cmp <- base::sum(!base::is.na(kh) & !base::is.na(ph)) + n_cmp <- sum(!is.na(kh) & !is.na(ph)) known_err[j] <- if (n_cmp == 0) NA_real_ else - base::round(base::sum(kh != ph, na.rm = TRUE) / n_cmp * 100, 2) + round(sum(kh != ph, na.rm = TRUE) / n_cmp * 100, 2) } } data.table::set(results_dt, j = "search_mode", value = search_mode) @@ -455,19 +451,19 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, "female_parent", known_id_vec[search_mode == "known_female"]) } - results_list <- base::vector("list", n_all) - for (j in base::seq_len(n_all)) { + results_list <- vector("list", n_all) + for (j in seq_len(n_all)) { prog_id <- progeny_ids[j] mode_j <- search_mode[j] if (mode_j == "pair") { pairs_j <- parent_pairs - pj <- base::match(j, pair_idx) + pj <- match(j, pair_idx) percent_mismatch <- pct_mismatch_mat[, pj] comparisons <- comparison_mat[, pj] } else { known_id <- known_id_vec[j] - if (!known_id %in% base::rownames(genos_mat)) { + if (!known_id %in% rownames(genos_mat)) { warning("Known parent '", known_id, "' of progeny '", prog_id, "' is not in the genotype file; this progeny cannot be tested.", call. = FALSE) @@ -479,8 +475,8 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, male_parent_candidates$id if (!allow_parent_selfing) cand <- cand[cand != known_id] if (exclude_self_match) cand <- cand[cand != prog_id] - cand <- cand[cand %in% base::rownames(genos_mat)] - if (base::length(cand) == 0) { + cand <- cand[cand %in% rownames(genos_mat)] + if (length(cand) == 0) { data.table::set(results_dt, j, "markers_tested", 0L) data.table::set(results_dt, j, "status", "low_markers") next @@ -490,12 +486,12 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, data.table::data.table(male_parent = cand, female_parent = known_id) pm <- genos_mat[pairs_j$male_parent, , drop = FALSE] pf <- genos_mat[pairs_j$female_parent, , drop = FALSE] - pg <- base::matrix(genos_mat[prog_id, ], nrow = base::nrow(pairs_j), - ncol = base::ncol(genos_mat), byrow = TRUE) - mm <- base::rowSums(.mend_mismatch(pm, pf, pg, ploidy), na.rm = TRUE) - comparisons <- base::rowSums(.mend_testable(pm, pf, pg, ploidy)) + pg <- matrix(genos_mat[prog_id, ], nrow = nrow(pairs_j), + ncol = ncol(genos_mat), byrow = TRUE) + mm <- rowSums(.mend_mismatch(pm, pf, pg, ploidy), na.rm = TRUE) + comparisons <- rowSums(.mend_testable(pm, pf, pg, ploidy)) percent_mismatch <- (mm / comparisons) * 100 - percent_mismatch[base::is.nan(percent_mismatch)] <- NA + percent_mismatch[is.nan(percent_mismatch)] <- NA } if (exclude_self_match) { @@ -504,9 +500,9 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, percent_mismatch[self_pair_idx] <- NA_real_ } - min_mismatch_val <- base::min(percent_mismatch, na.rm = TRUE) + min_mismatch_val <- min(percent_mismatch, na.rm = TRUE) - if (base::is.infinite(min_mismatch_val)) { + if (is.infinite(min_mismatch_val)) { data.table::set(results_dt, j, "markers_tested", 0L) data.table::set(results_dt, j, "status", "low_markers") next @@ -514,64 +510,64 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, if (top_n > 1) { # Rank all candidate pairs by error, then by markers tested; keep top_n - valid_idx <- base::which(!base::is.na(percent_mismatch)) - ranked_idx <- valid_idx[base::order(percent_mismatch[valid_idx], + valid_idx <- which(!is.na(percent_mismatch)) + ranked_idx <- valid_idx[order(percent_mismatch[valid_idx], -comparisons[valid_idx])] - best_indices <- ranked_idx[base::seq_len(base::min(top_n, base::length(ranked_idx)))] + best_indices <- ranked_idx[seq_len(min(top_n, length(ranked_idx)))] } else { - best_indices <- base::which(percent_mismatch == min_mismatch_val) - if (base::length(best_indices) > 1) { + best_indices <- which(percent_mismatch == min_mismatch_val) + if (length(best_indices) > 1) { best_markers_per_pair <- comparisons[best_indices] - max_markers <- base::max(best_markers_per_pair) + max_markers <- max(best_markers_per_pair) best_indices <- best_indices[best_markers_per_pair == max_markers] } } best_pairs <- pairs_j[best_indices] best_markers <- comparisons[best_indices[1]] - best_error <- base::round(min_mismatch_val, 2) + best_error <- round(min_mismatch_val, 2) a_status <- assign_status(best_markers, best_error) - if (top_n == 1 && !show_ties && base::nrow(best_pairs) > 1) { - warning("Progeny '", prog_id, "' has ", base::nrow(best_pairs), + if (top_n == 1 && !show_ties && nrow(best_pairs) > 1) { + warning("Progeny '", prog_id, "' has ", nrow(best_pairs), " tied best pairs. Only one is reported as show_ties=FALSE.", call. = FALSE) } report_all <- show_ties || top_n > 1 - num_to_report <- base::min(base::nrow(best_pairs), - if (report_all) base::nrow(best_pairs) else 1) + num_to_report <- min(nrow(best_pairs), + if (report_all) nrow(best_pairs) else 1) data.table::set(results_dt, j, "male_parent", best_pairs$male_parent[1]) data.table::set(results_dt, j, "female_parent", best_pairs$female_parent[1]) - data.table::set(results_dt, j, "mendelian_error_pct", base::sprintf("%.2f", min_mismatch_val)) - data.table::set(results_dt, j, "markers_tested", base::as.integer(best_markers)) + data.table::set(results_dt, j, "mendelian_error_pct", sprintf("%.2f", min_mismatch_val)) + data.table::set(results_dt, j, "markers_tested", as.integer(best_markers)) data.table::set(results_dt, j, "status", a_status) data.table::set(results_dt, j, "tied_candidates", - base::as.integer(base::sum(percent_mismatch == min_mismatch_val, + as.integer(sum(percent_mismatch == min_mismatch_val, na.rm = TRUE))) if (report_all && num_to_report > 1) { - tie_row <- base::list(id = prog_id) - for (k in base::seq(2, num_to_report)) { - tie_row[[base::paste0("male_parent_", k)]] <- best_pairs$male_parent[k] - tie_row[[base::paste0("female_parent_", k)]] <- best_pairs$female_parent[k] - tie_row[[base::paste0("mendelian_error_pct_", k)]] <- - if (top_n > 1) base::round(percent_mismatch[best_indices[k]], 2) else min_mismatch_val - tie_row[[base::paste0("markers_tested_", k)]] <- comparisons[best_indices[k]] + tie_row <- list(id = prog_id) + for (k in seq(2, num_to_report)) { + tie_row[[paste0("male_parent_", k)]] <- best_pairs$male_parent[k] + tie_row[[paste0("female_parent_", k)]] <- best_pairs$female_parent[k] + tie_row[[paste0("mendelian_error_pct_", k)]] <- + if (top_n > 1) round(percent_mismatch[best_indices[k]], 2) else min_mismatch_val + tie_row[[paste0("markers_tested_", k)]] <- comparisons[best_indices[k]] } results_list[[j]] <- data.table::as.data.table(tie_row) } } tie_rows <- data.table::rbindlist( - base::Filter(Negate(base::is.null), results_list), + Filter(Negate(is.null), results_list), fill = TRUE, use.names = TRUE ) - if (base::nrow(tie_rows) > 0) { + if (nrow(tie_rows) > 0) { final_df <- merge(results_dt, tie_rows, by = "id", all.x = TRUE) - for (col in base::names(final_df)) + for (col in names(final_df)) data.table::set(final_df, which(final_df[[col]] == ""), col, NA_character_) } else { final_df <- results_dt @@ -589,40 +585,39 @@ find_parentage <- function(genotypes_file, parents_file, progeny_file, #### Verbose output #### if (verbose) { - total_progeny <- base::nrow(final_df) - base::cat("\n=== Parentage Assignment Report ===\n") - base::cat("\nTotal progeny evaluated:", total_progeny, "\n") - base::cat("Method:", method_reported, " | ", - "Error threshold:", error_threshold, "% | ", - "Minimum markers:", min_markers, "\n") + total_progeny <- nrow(final_df) + cat("\n=== Parentage Assignment Report ===\n") + cat(sprintf("\nTotal progeny evaluated: %d\n", total_progeny)) + cat(sprintf("Method: %s | Error threshold: %s%% | Minimum markers: %s\n", + method_reported, error_threshold, min_markers)) - n_pass <- base::nrow(output_list$pass) + n_pass <- nrow(output_list$pass) if (n_pass > 0) { - base::cat(base::sprintf("\n%d progeny passed (%.1f%%).\n", + cat(sprintf("\n%d progeny passed (%.1f%%).\n", n_pass, (n_pass / total_progeny) * 100)) } else { - base::cat("\nNo progeny passed.\n") + cat("\nNo progeny passed.\n") } - n_high <- base::nrow(output_list$high_error) + n_high <- nrow(output_list$high_error) if (n_high > 0) { - base::cat(base::sprintf("\n%d progeny flagged high_error (%.1f%%):\n", + cat(sprintf("\n%d progeny flagged high_error (%.1f%%):\n", n_high, (n_high / total_progeny) * 100)) - base::print(output_list$high_error) + print(output_list$high_error) } else { - base::cat("\nNo progeny flagged for high error.\n") + cat("\nNo progeny flagged for high error.\n") } - n_low <- base::nrow(output_list$low_markers) + n_low <- nrow(output_list$low_markers) if (n_low > 0) { - base::cat(base::sprintf("\n%d progeny flagged low_markers (%.1f%%):\n", + cat(sprintf("\n%d progeny flagged low_markers (%.1f%%):\n", n_low, (n_low / total_progeny) * 100)) - base::print(output_list$low_markers) + print(output_list$low_markers) } else { - base::cat("\nNo progeny flagged for low marker count.\n") + cat("\nNo progeny flagged for low marker count.\n") } - base::cat("\nFull results are included in the returned list as $full_results.\n") + cat("\nFull results are included in the returned list as $full_results.\n") } #### Plot Results #### @@ -632,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, " | ", @@ -673,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", @@ -682,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 6027216..ca73095 100644 --- a/R/utils.R +++ b/R/utils.R @@ -25,8 +25,8 @@ utils::globalVariables(c( #' @param ploidy even integer ploidy level. #' @return Object matching \code{g} holding the gamete dosage bound. #' @noRd -.gamete_lo <- function(g, ploidy) base::pmax(0L, g - ploidy / 2L) -.gamete_hi <- function(g, ploidy) base::pmin(ploidy / 2L, g) +.gamete_lo <- function(g, ploidy) pmax(0L, g - ploidy / 2L) +.gamete_hi <- function(g, ploidy) pmin(ploidy / 2L, g) #' Flag impossible offspring dosages #' @@ -79,10 +79,10 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @noRd .mend_testable <- function(male, female, offspring, ploidy) { if (ploidy %% 2 == 0) - return(!base::is.na(male) & !base::is.na(female) & !base::is.na(offspring)) - o_hom <- !base::is.na(offspring) & (offspring == 0 | offspring == ploidy) - m_hom <- !base::is.na(male) & (male == 0 | male == ploidy) - f_hom <- !base::is.na(female) & (female == 0 | female == ploidy) + return(!is.na(male) & !is.na(female) & !is.na(offspring)) + o_hom <- !is.na(offspring) & (offspring == 0 | offspring == ploidy) + m_hom <- !is.na(male) & (male == 0 | male == ploidy) + f_hom <- !is.na(female) & (female == 0 | female == ploidy) o_hom & (m_hom | f_hom) } @@ -94,10 +94,10 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @return Invisibly TRUE if valid; otherwise an error is thrown. #' @noRd .check_ploidy <- function(ploidy) { - if (base::length(ploidy) != 1 || !base::is.numeric(ploidy) || - base::is.na(ploidy) || ploidy < 2 || ploidy != base::round(ploidy)) - base::stop("ploidy must be an integer >= 2.") - base::invisible(TRUE) + if (length(ploidy) != 1 || !is.numeric(ploidy) || + is.na(ploidy) || ploidy < 2 || ploidy != round(ploidy)) + stop("ploidy must be an integer >= 2.") + invisible(TRUE) } #### Column name helpers #### @@ -112,7 +112,7 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @return Normalized character vector. #' @noRd .normalize_name <- function(x) { - base::gsub("[ .-]+", "_", base::tolower(base::trimws(x))) + gsub("[ .-]+", "_", tolower(trimws(x))) } #' Rename columns to their standard names, ignoring case @@ -127,23 +127,23 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @return \code{x} with matching columns renamed. #' @noRd .standardize_names <- function(x, cols) { - nm <- base::names(x) + nm <- names(x) norm <- .normalize_name(nm) - old <- new <- base::character(0) + old <- new <- character(0) for (col in cols) { if (col %in% nm) next - hit <- base::which(norm == col & !(nm %in% old)) - if (base::length(hit) >= 1) { + hit <- which(norm == col & !(nm %in% old)) + if (length(hit) >= 1) { old <- c(old, nm[hit[1]]) new <- c(new, col) } } - if (base::length(old) == 0) return(x) + if (length(old) == 0) return(x) if (data.table::is.data.table(x)) { x <- data.table::copy(x) data.table::setnames(x, old, new) } else { - base::names(x)[base::match(old, base::names(x))] <- new + names(x)[match(old, names(x))] <- new } x } @@ -160,15 +160,15 @@ mendelian_error <- function(male, female, offspring, ploidy) { #' @return One of "vcf", "plink", "text", "data.frame", or "matrix". #' @noRd .genotype_source <- function(x) { - if (base::inherits(x, "vcfR")) return("vcf") - if (base::is.matrix(x)) return("matrix") - if (base::is.data.frame(x)) return("data.frame") - if (base::is.character(x) && base::length(x) == 1) { - if (base::grepl("\\.vcf(\\.gz)?$", x, ignore.case = TRUE)) return("vcf") - if (base::grepl("\\.ped$", x, ignore.case = TRUE)) return("plink") + if (inherits(x, "vcfR")) return("vcf") + if (is.matrix(x)) return("matrix") + if (is.data.frame(x)) return("data.frame") + if (is.character(x) && length(x) == 1) { + if (grepl("\\.vcf(\\.gz)?$", x, ignore.case = TRUE)) return("vcf") + if (grepl("\\.ped$", x, ignore.case = TRUE)) return("plink") return("text") } - base::stop("Genotypes must be a file path (TXT/TSV/CSV, VCF or PLINK .ped), ", + stop("Genotypes must be a file path (TXT/TSV/CSV, VCF or PLINK .ped), ", "a vcfR object, a data.frame / data.table, or a matrix.") } @@ -196,7 +196,7 @@ mendelian_error <- function(male, female, offspring, ploidy) { format = c("data.frame", "matrix"), verbose = FALSE, counted_allele = NULL) { - format <- base::match.arg(format) + format <- match.arg(format) src <- .genotype_source(x) # VCF: convert straight to the requested layout @@ -205,8 +205,8 @@ mendelian_error <- function(male, female, offspring, ploidy) { # PLINK .ped: diploid only; a .map next to the file supplies marker names if (src == "plink") { - if (!base::is.null(ploidy) && ploidy != 2) - base::stop("PLINK .ped files are diploid; use ploidy = 2.") + if (!is.null(ploidy) && ploidy != 2) + stop("PLINK .ped files are diploid; use ploidy = 2.") return(ped_to_dosage(x, format = format, counted_allele = counted_allele, verbose = verbose)) } @@ -214,39 +214,39 @@ mendelian_error <- function(male, female, offspring, ploidy) { # Text file: read to a table, then fall through to the in-memory handling from_text <- src == "text" if (from_text) { - if (!base::file.exists(x)) - base::stop("Genotype file not found: ", x) + if (!file.exists(x)) + stop("Genotype file not found: ", x) x <- data.table::fread(x, sep = "auto", data.table = FALSE) src <- "data.frame" } # ID column: "id" in any case (id, ID, Id, ...); exact "id" preferred id_col <- if (src == "data.frame") { - nm <- base::names(x) + nm <- names(x) hit <- if ("id" %in% nm) "id" else nm[.normalize_name(nm) == "id"] - if (base::length(hit) >= 1) hit[1] else NA_character_ + if (length(hit) >= 1) hit[1] else NA_character_ } else { NA_character_ } # Text files have no row names, so the first column holds the IDs - if (from_text && format == "matrix" && base::is.na(id_col)) - id_col <- base::names(x)[1] + if (from_text && format == "matrix" && is.na(id_col)) + id_col <- names(x)[1] if (format == "data.frame") { if (src == "matrix") return(data.table::as.data.table(x, keep.rownames = "id")) out <- data.table::as.data.table(x) - if (!base::is.na(id_col) && id_col != "id") + if (!is.na(id_col) && id_col != "id") data.table::setnames(out, id_col, "id") return(out) } # format == "matrix": individuals in named rows, markers in columns - if (src == "matrix" || base::is.na(id_col)) return(x) - x <- base::as.data.frame(x) # data.table subsetting differs from data.frame - ids <- base::as.character(x[[id_col]]) - x <- base::as.matrix(x[, base::names(x) != id_col, drop = FALSE]) - base::rownames(x) <- ids + if (src == "matrix" || is.na(id_col)) return(x) + x <- as.data.frame(x) # data.table subsetting differs from data.frame + ids <- as.character(x[[id_col]]) + x <- as.matrix(x[, names(x) != id_col, drop = FALSE]) + rownames(x) <- ids x } @@ -323,7 +323,7 @@ QPsolve <- function(Y, X) { if (ploidy %% 2 == 0) { h <- ploidy / 2 parent_bad <- function(g) { - !(offspring >= base::pmax(0, g - h) & offspring <= base::pmin(g, h) + h) + !(offspring >= pmax(0, g - h) & offspring <= pmin(g, h) + h) } male_bad <- parent_bad(male) female_bad <- parent_bad(female) @@ -334,11 +334,11 @@ QPsolve <- function(Y, X) { } # missing calls are named after every individual without a call, e.g. # missing_male, missing_male_female, missing_progeny_female - na_p <- base::is.na(offspring); na_m <- base::is.na(male); na_f <- base::is.na(female) - missing_label <- base::paste0("missing", - base::ifelse(na_p, "_progeny", ""), - base::ifelse(na_m, "_male", ""), - base::ifelse(na_f, "_female", "")) + na_p <- is.na(offspring); na_m <- is.na(male); na_f <- is.na(female) + missing_label <- paste0("missing", + ifelse(na_p, "_progeny", ""), + ifelse(na_m, "_male", ""), + ifelse(na_f, "_female", "")) dplyr::case_when( na_p | na_m | na_f ~ missing_label, !tt ~ "untestable", @@ -349,3 +349,25 @@ QPsolve <- function(Y, X) { TRUE ~ "mismatch_combination" ) } + + +# Read founder ids from a file path, a character vector of ids, or a +# data.frame / data.table (first column). Returns a character vector. +.read_founders <- function(x) { + if (is.null(x)) return(character(0)) + if (is.data.frame(x)) return(unique(as.character(x[[1]]))) + if (is.character(x) && length(x) == 1 && file.exists(x)) { + out <- tryCatch( + data.table::fread(x, header = FALSE, colClasses = "character"), + error = function(e) + stop("Could not read founders list. Ensure it is a plain text or CSV/TSV file.") + ) + return(unique(out[[1]])) + } + # 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 ea723fb..56805ad 100644 --- a/R/validate_pedigree.R +++ b/R/validate_pedigree.R @@ -15,8 +15,8 @@ #' (diploid only), converted with [ped_to_dosage()]; or a data.frame / #' data.table with an id column followed by marker columns coded as #' allele-B dosage (0, 1, ..., ploidy). -#' @param founders_file Character, optional. Path to a one-column file listing -#' founder IDs. Founders with both parents coded as 0 are left unchanged. +#' @param founders_file Optional. Founder IDs as a path to a one-column file, a +#' character vector of ids, or a data.frame (first column). Founders with both parents coded as 0 are left unchanged. #' Defaults to NULL. #' @param trio_error_threshold Numeric. Maximum Mendelian error percentage to #' classify a trio as pass (default: 5.0). Must be between 0 and 100. @@ -178,11 +178,11 @@ validate_pedigree <- function(pedigree_file, genotypes_file, # Column names are matched ignoring case (e.g. ID, Male_Parent, FEMALE_PARENT) pedigree <- .standardize_names(pedigree, c("id", "male_parent", "female_parent")) required_ped_cols <- c("id", "male_parent", "female_parent") - missing_cols <- base::setdiff(required_ped_cols, base::names(pedigree)) - if (base::length(missing_cols) > 0) + missing_cols <- setdiff(required_ped_cols, names(pedigree)) + if (length(missing_cols) > 0) stop("Pedigree file missing required columns: ", - base::paste(missing_cols, collapse = ", ")) - if (!"id" %in% base::names(genos)) + paste(missing_cols, collapse = ", ")) + if (!"id" %in% names(genos)) stop("Genotypes file must have an 'id' column") pedigree[, male_parent := as.character(male_parent)] @@ -190,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) @@ -217,20 +208,20 @@ validate_pedigree <- function(pedigree_file, genotypes_file, no_geno_rows <- pedigree[!(id %in% valid_ids) | (!(male_parent %in% valid_ids) & male_parent != "0") | (!(female_parent %in% valid_ids) & female_parent != "0")] - if (base::nrow(no_geno_rows) > 0 && verbose) - base::cat("Found", base::nrow(no_geno_rows), - "trios with missing genotype data; flagged as no_genotype_data.\n") + if (nrow(no_geno_rows) > 0 && verbose) + cat(sprintf("Found %d trios with missing genotype data; flagged as no_genotype_data.\n", + nrow(no_geno_rows))) pedigree <- has_geno - if (base::nrow(pedigree) == 0) + if (nrow(pedigree) == 0) stop("No valid trios remain after filtering for genotype availability.") #### Per-marker mismatch counters (filled in the trio loop) #### - marker_names <- base::colnames(genos_mat) - marker_tested_n <- base::numeric(base::length(marker_names)) - marker_mismatch_n <- base::numeric(base::length(marker_names)) + marker_names <- colnames(genos_mat) + marker_tested_n <- numeric(length(marker_names)) + marker_mismatch_n <- numeric(length(marker_names)) #### Main trio evaluation loop #### - results_list <- base::lapply(base::seq_len(base::nrow(pedigree)), function(i) { + results_list <- lapply(seq_len(nrow(pedigree)), function(i) { prog_id <- pedigree$id[i] male_parent_id <- pedigree$male_parent[i] female_parent_id <- pedigree$female_parent[i] @@ -263,11 +254,11 @@ validate_pedigree <- function(pedigree_file, genotypes_file, progeny_vec, ploidy) tt_vec <- .mend_testable(male_parent_vec, female_parent_vec, progeny_vec, ploidy) - mismatches <- base::sum(mm_vec, na.rm = TRUE) - markers_tested <- base::sum(tt_vec) - marker_tested_n <<- marker_tested_n + base::as.numeric(tt_vec) + mismatches <- sum(mm_vec, na.rm = TRUE) + markers_tested <- sum(tt_vec) + marker_tested_n <<- marker_tested_n + as.numeric(tt_vec) marker_mismatch_n <<- marker_mismatch_n + - base::as.numeric((mm_vec %in% TRUE) & tt_vec) + as.numeric((mm_vec %in% TRUE) & tt_vec) if (markers_tested == 0) { status <- "no_data" correction_decision <- "none" @@ -285,15 +276,15 @@ validate_pedigree <- function(pedigree_file, genotypes_file, progeny_hom <- genos_hom_mat[prog_id, ] male_parent_hom <- genos_hom_mat[male_parent_id, ] female_parent_hom <- genos_hom_mat[female_parent_id, ] - male_comparisons <- base::sum(!base::is.na(male_parent_hom) & - !base::is.na(progeny_hom)) + male_comparisons <- sum(!is.na(male_parent_hom) & + !is.na(progeny_hom)) male_parent_error_pct <- if (male_comparisons == 0) NA_real_ else - base::round((base::sum(male_parent_hom != progeny_hom, na.rm = TRUE) / + round((sum(male_parent_hom != progeny_hom, na.rm = TRUE) / male_comparisons) * 100, 2) - female_comparisons <- base::sum(!base::is.na(female_parent_hom) & - !base::is.na(progeny_hom)) + female_comparisons <- sum(!is.na(female_parent_hom) & + !is.na(progeny_hom)) female_parent_error_pct <- if (female_comparisons == 0) NA_real_ else - base::round((base::sum(female_parent_hom != progeny_hom, na.rm = TRUE) / + round((sum(female_parent_hom != progeny_hom, na.rm = TRUE) / female_comparisons) * 100, 2) male_acceptable <- !is.na(male_parent_error_pct) && male_parent_error_pct <= single_parent_error_threshold @@ -303,7 +294,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, # Neither parent is excluded by the homozygous check. If the trio # still fails (e.g. a recorded self whose offspring carries alleles # absent from the parent), the pedigree is wrong but the faulty - # parent cannot be identified: flag it instead of keeping it. + # parent cannot be identified, so the trio is flagged and the parents are not kept. correction_decision <- if (error_pct > trio_error_threshold) "unresolved" else "keep_both" } else if (male_acceptable && !female_acceptable) { correction_decision <- "remove_female_parent" @@ -323,7 +314,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, id = prog_id, orig_male_parent = male_parent_id, orig_female_parent = female_parent_id, - trio_mendelian_error_pct = base::round(error_pct, 2), + trio_mendelian_error_pct = round(error_pct, 2), trio_markers_tested = markers_tested, status = status, recommended_correction = correction_decision, @@ -335,7 +326,7 @@ validate_pedigree <- function(pedigree_file, genotypes_file, final_df <- data.table::rbindlist(results_list) #### Append no_genotype_data rows #### - if (base::nrow(no_geno_rows) > 0) { + if (nrow(no_geno_rows) > 0) { no_geno_df <- data.table::data.table( id = no_geno_rows$id, orig_male_parent = no_geno_rows$male_parent, @@ -352,12 +343,12 @@ validate_pedigree <- function(pedigree_file, genotypes_file, #### Build corrected pedigree #### corrected_pedigree <- data.table::copy(original_pedigree) - for (i in base::seq_len(base::nrow(final_df))) { + for (i in seq_len(nrow(final_df))) { prog_id <- final_df$id[i] decision <- final_df$recommended_correction[i] - row_idx <- base::which(corrected_pedigree$id == prog_id) + row_idx <- which(corrected_pedigree$id == prog_id) # low_markers_* decisions apply the same removal as the plain decision - base_decision <- base::sub("^low_markers_", "", decision) + base_decision <- sub("^low_markers_", "", decision) if (base_decision %in% c("remove_male_parent", "remove_both")) data.table::set(corrected_pedigree, row_idx, "male_parent", "0") if (base_decision %in% c("remove_female_parent", "remove_both")) @@ -366,24 +357,23 @@ validate_pedigree <- function(pedigree_file, genotypes_file, #### Summary output #### if (verbose) { - total_trios <- base::nrow(final_df) - status_counts <- base::table(final_df$status) - base::cat("\n--- Trio Validation Summary ---\n") - base::cat("Total trios in pedigree:", total_trios, "\n") - for (s in base::names(status_counts)) - base::cat(base::sprintf("%-28s: %d (%.1f%%)\n", s, - status_counts[s], - (status_counts[s] / total_trios) * 100)) - base::cat("Error threshold:", trio_error_threshold, "%\n") - base::cat("Homozygous threshold:", single_parent_error_threshold, "%\n") - base::cat("Minimum markers required:", min_markers, "\n\n") - corrections <- base::table(final_df$recommended_correction) - base::cat("Correction summary:\n") - for (decision in base::names(corrections)) + total_trios <- nrow(final_df) + status_counts <- table(final_df$status) + cat("\n--- Trio Validation Summary ---\n") + cat(sprintf("Total trios in pedigree: %d\n", total_trios)) + for (s in names(status_counts)) + cat(sprintf("%-28s: %d (%.1f%%)\n", s, status_counts[s], + (status_counts[s] / total_trios) * 100)) + cat(sprintf("Error threshold: %s%%\n", trio_error_threshold)) + cat(sprintf("Homozygous threshold: %s%%\n", single_parent_error_threshold)) + cat(sprintf("Minimum markers required: %s\n\n", min_markers)) + corrections <- table(final_df$recommended_correction) + cat("Correction summary:\n") + for (decision in names(corrections)) if (decision != "none") - base::cat(" ", decision, ":", corrections[decision], "\n") - base::cat("\n") - base::print(final_df) + cat(sprintf(" %s: %d\n", decision, corrections[decision])) + cat("\n") + print(final_df) } #### Plot results #### @@ -392,20 +382,21 @@ validate_pedigree <- function(pedigree_file, genotypes_file, if (!requireNamespace("ggplot2", quietly = TRUE)) { warning("ggplot2 is required for plot_results = TRUE. Please install it.", call. = FALSE) } else { - plot_df <- final_df[!is.na(final_df$trio_mendelian_error_pct)] # Colour by the trio error first, so no trio above the threshold is shown as pass - plot_df$plot_status <- dplyr::case_when( - plot_df$trio_mendelian_error_pct <= trio_error_threshold ~ "pass", - plot_df$recommended_correction %in% c("remove_male_parent", - "remove_female_parent", - "low_markers_remove_male_parent", - "low_markers_remove_female_parent") ~ "fail_one_parent", - plot_df$recommended_correction %in% c("remove_both", - "low_markers_remove_both") ~ "fail_both_parents", - plot_df$recommended_correction %in% c("unresolved", - "low_markers_unresolved") ~ "fail_unresolved", - TRUE ~ "other" - ) + plot_df <- final_df %>% + dplyr::filter(!is.na(trio_mendelian_error_pct)) %>% + dplyr::mutate(plot_status = dplyr::case_when( + trio_mendelian_error_pct <= trio_error_threshold ~ "pass", + recommended_correction %in% c("remove_male_parent", + "remove_female_parent", + "low_markers_remove_male_parent", + "low_markers_remove_female_parent") ~ "fail_one_parent", + recommended_correction %in% c("remove_both", + "low_markers_remove_both") ~ "fail_both_parents", + recommended_correction %in% c("unresolved", + "low_markers_unresolved") ~ "fail_unresolved", + TRUE ~ "other" + )) n_total <- nrow(plot_df) n_fail <- sum(plot_df$trio_mendelian_error_pct > trio_error_threshold) n_pass <- sum(plot_df$trio_mendelian_error_pct <= trio_error_threshold) @@ -453,29 +444,30 @@ validate_pedigree <- function(pedigree_file, genotypes_file, marker_trio_out <- NULL if (marker_trio_table) { full_trios <- pedigree[male_parent != "0" & female_parent != "0"] - status_list <- base::lapply(base::seq_len(base::nrow(full_trios)), function(i) + status_list <- lapply(seq_len(nrow(full_trios)), function(i) .marker_trio_status(genos_mat[full_trios$male_parent[i], ], genos_mat[full_trios$female_parent[i], ], genos_mat[full_trios$id[i], ], ploidy)) - base::names(status_list) <- base::make.unique( - base::paste(full_trios$id, full_trios$male_parent, full_trios$female_parent, sep = "_")) + names(status_list) <- make.unique( + paste(full_trios$id, full_trios$male_parent, full_trios$female_parent, sep = "_")) marker_trio_out <- data.table::as.data.table( - base::c(base::list(marker = marker_names), status_list)) + c(list(marker = marker_names), status_list)) } #### Per-marker mismatch summary #### marker_summary <- data.table::data.table( marker = marker_names, - trios_tested = base::as.integer(marker_tested_n), - trios_mismatch = base::as.integer(marker_mismatch_n) - ) - marker_summary[, mismatch_pct := base::ifelse( - trios_tested > 0, base::round(trios_mismatch / trios_tested * 100, 2), NA_real_)] - marker_summary <- marker_summary[base::order(-mismatch_pct, -trios_mismatch, marker)] + trios_tested = as.integer(marker_tested_n), + trios_mismatch = as.integer(marker_mismatch_n) + ) %>% + dplyr::mutate(mismatch_pct = ifelse(trios_tested > 0, + round(trios_mismatch / trios_tested * 100, 2), + NA_real_)) %>% + dplyr::arrange(dplyr::desc(mismatch_pct), dplyr::desc(trios_mismatch), marker) #### Compile and return named list #### - output_list <- base::list( + output_list <- list( pass = final_df[status == "pass"], fail = final_df[status == "fail"], low_markers = final_df[status == "low_markers"], @@ -490,5 +482,5 @@ validate_pedigree <- function(pedigree_file, genotypes_file, marker_trio_table = marker_trio_out, plot = p ) - return(base::invisible(output_list)) + return(invisible(output_list)) } diff --git a/R/vcf_to_dosage.R b/R/vcf_to_dosage.R index 7d753cf..149f945 100644 --- a/R/vcf_to_dosage.R +++ b/R/vcf_to_dosage.R @@ -107,27 +107,27 @@ vcf_to_dosage <- function(vcf, ...) { #### Input Validation #### - format <- base::match.arg(format) - marker_names <- base::match.arg(marker_names) + format <- match.arg(format) + marker_names <- match.arg(marker_names) - if (!base::requireNamespace("vcfR", quietly = TRUE)) + if (!requireNamespace("vcfR", quietly = TRUE)) stop("Package 'vcfR' is required to read VCF files. Install it with install.packages(\"vcfR\").") - if (!base::is.null(ploidy)) .check_ploidy(ploidy) - if (!base::is.null(dosage_field) && base::is.null(ploidy)) + if (!is.null(ploidy)) .check_ploidy(ploidy) + if (!is.null(dosage_field) && is.null(ploidy)) stop("ploidy must be supplied when dosage_field is used; it cannot be inferred from a dosage field.") for (arg in c("min_marker_call_rate", "min_sample_call_rate")) { - val <- base::get(arg) - if (base::length(val) != 1 || !base::is.numeric(val) || base::is.na(val) || + val <- get(arg) + if (length(val) != 1 || !is.numeric(val) || is.na(val) || val < 0 || val > 1) stop(arg, " must be a single number between 0 and 1.") } #### Read the VCF #### - if (base::is.character(vcf) && base::length(vcf) == 1) { - if (!base::file.exists(vcf)) + if (is.character(vcf) && length(vcf) == 1) { + if (!file.exists(vcf)) stop("VCF file not found: ", vcf) vcf_obj <- vcfR::read.vcfR(vcf, verbose = verbose, ...) - } else if (base::inherits(vcf, "vcfR")) { + } else if (inherits(vcf, "vcfR")) { vcf_obj <- vcf } else { stop("vcf must be a path to a .vcf / .vcf.gz file or a vcfR object.") @@ -136,23 +136,23 @@ vcf_to_dosage <- function(vcf, # getFIX() drops to a named vector when the VCF has one variant; keep a matrix get_fix <- function(obj) { f <- vcfR::getFIX(obj) - if (!base::is.null(f) && base::is.null(base::dim(f))) - f <- base::matrix(f, nrow = 1, dimnames = base::list(NULL, base::names(f))) + if (!is.null(f) && is.null(dim(f))) + f <- matrix(f, nrow = 1, dimnames = list(NULL, names(f))) f } fix <- get_fix(vcf_obj) - if (base::is.null(fix) || base::nrow(fix) == 0) + if (is.null(fix) || nrow(fix) == 0) stop("The VCF contains no variants.") #### Drop multiallelic markers #### - n_markers_in <- base::nrow(fix) + n_markers_in <- nrow(fix) if (biallelic_only) { - keep <- !base::grepl(",", fix[, "ALT"], fixed = TRUE) & !base::is.na(fix[, "ALT"]) - if (base::sum(keep) == 0) + keep <- !grepl(",", fix[, "ALT"], fixed = TRUE) & !is.na(fix[, "ALT"]) + if (sum(keep) == 0) stop("No biallelic markers remain. Set biallelic_only = FALSE to keep multiallelic markers.") - if (base::any(!keep)) { - warning(base::sum(!keep), " multiallelic marker(s) were dropped. ", + if (any(!keep)) { + warning(sum(!keep), " multiallelic marker(s) were dropped. ", "Set biallelic_only = FALSE to keep them.", call. = FALSE) vcf_obj <- vcf_obj[keep, ] fix <- get_fix(vcf_obj) @@ -161,9 +161,9 @@ vcf_to_dosage <- function(vcf, #### Marker names #### ids <- fix[, "ID"] - positions <- base::paste(fix[, "CHROM"], fix[, "POS"], sep = "_") - id_usable <- !base::any(base::is.na(ids) | ids == "." | ids == "") && - !base::any(base::duplicated(ids)) + positions <- paste(fix[, "CHROM"], fix[, "POS"], sep = "_") + id_usable <- !any(is.na(ids) | ids == "." | ids == "") && + !any(duplicated(ids)) if (marker_names == "id") { if (!id_usable) @@ -174,21 +174,21 @@ vcf_to_dosage <- function(vcf, } else { snp_names <- if (id_usable) ids else positions } - if (base::any(base::duplicated(snp_names))) { + if (any(duplicated(snp_names))) { warning("Duplicated marker names were made unique with a numeric suffix.", call. = FALSE) - snp_names <- base::make.unique(snp_names, sep = "_") + snp_names <- make.unique(snp_names, sep = "_") } #### Extract genotype calls #### - if (!base::is.null(dosage_field)) { + if (!is.null(dosage_field)) { # Dosages are already numeric in the requested FORMAT field dose <- vcfR::extract.gt(vcf_obj, element = dosage_field, as.numeric = TRUE, IDtoRowNames = FALSE) - if (base::all(base::is.na(dose))) + if (all(is.na(dose))) stop("FORMAT field '", dosage_field, "' is absent or empty in this VCF.") - obs_max <- base::suppressWarnings(base::max(dose, na.rm = TRUE)) - if (base::is.finite(obs_max) && obs_max > ploidy) + obs_max <- suppressWarnings(max(dose, na.rm = TRUE)) + if (is.finite(obs_max) && obs_max > ploidy) warning("Dosages above the requested ploidy (", ploidy, ") were found in field '", dosage_field, "'.", call. = FALSE) n_wrong_ploidy <- 0L @@ -196,47 +196,47 @@ vcf_to_dosage <- function(vcf, } else { gt <- vcfR::extract.gt(vcf_obj, element = "GT", IDtoRowNames = FALSE) - if (base::is.null(gt) || base::all(base::is.na(gt))) + if (is.null(gt) || all(is.na(gt))) stop("No GT field found in this VCF.") # Genotype strings repeat heavily, so the parsing is done once per unique # call and mapped back to the full matrix. This keeps large files cheap. - calls <- base::unique(base::as.vector(gt)) - idx <- base::match(gt, calls) - n_obs <- base::tabulate(idx, nbins = base::length(calls)) + calls <- unique(as.vector(gt)) + idx <- match(gt, calls) + n_obs <- tabulate(idx, nbins = length(calls)) - alleles <- base::strsplit(calls, "[/|]") # phased or unphased - n_slots <- base::lengths(alleles) - is_miss <- base::is.na(calls) | n_slots == 0 | - base::vapply(alleles, - function(a) base::any(base::is.na(a) | a == "."), - base::logical(1)) + alleles <- strsplit(calls, "[/|]") # phased or unphased + n_slots <- lengths(alleles) + is_miss <- is.na(calls) | n_slots == 0 | + vapply(alleles, + function(a) any(is.na(a) | a == "."), + logical(1)) # dosage = allele slots that are not the reference allele (0) - dose_u <- base::rep(NA_real_, base::length(calls)) - dose_u[!is_miss] <- base::vapply(alleles[!is_miss], - function(a) base::sum(a != "0"), - base::numeric(1)) + dose_u <- rep(NA_real_, length(calls)) + dose_u[!is_miss] <- vapply(alleles[!is_miss], + function(a) sum(a != "0"), + numeric(1)) #### Ploidy #### obs_slots <- n_slots[!is_miss] obs_count <- n_obs[!is_miss] - if (base::length(obs_slots) == 0) + if (length(obs_slots) == 0) stop("Every genotype call is missing; nothing to convert.") - slot_totals <- base::tapply(obs_count, obs_slots, base::sum) - obs_ploidy <- base::as.integer(base::names(slot_totals)[base::which.max(slot_totals)]) + slot_totals <- tapply(obs_count, obs_slots, sum) + obs_ploidy <- as.integer(names(slot_totals)[which.max(slot_totals)]) ploidy_mismatch <- FALSE - if (base::is.null(ploidy)) { + if (is.null(ploidy)) { ploidy <- obs_ploidy if (verbose) - base::cat("Inferred ploidy from the GT field:", ploidy, "\n") + cat(sprintf("Inferred ploidy from the GT field: %s\n", ploidy)) } else { ploidy_mismatch <- obs_ploidy != ploidy } off_ploidy <- !is_miss & n_slots != ploidy - n_wrong_ploidy <- base::sum(n_obs[off_ploidy]) + n_wrong_ploidy <- sum(n_obs[off_ploidy]) # One warning per situation: a contradicting ploidy, or a few stray calls if (ploidy_mismatch) { @@ -250,58 +250,56 @@ vcf_to_dosage <- function(vcf, } dose_u[is_miss | off_ploidy] <- NA_real_ - dose <- base::matrix(dose_u[idx], - nrow = base::nrow(gt), ncol = base::ncol(gt), - dimnames = base::list(NULL, base::colnames(gt))) - base::storage.mode(dose) <- "integer" # dosages are whole numbers + dose <- matrix(dose_u[idx], + nrow = nrow(gt), ncol = ncol(gt), + dimnames = list(NULL, colnames(gt))) + storage.mode(dose) <- "integer" # dosages are whole numbers } #### Reference instead of alternate dosage #### - ploidy <- base::as.integer(ploidy) + ploidy <- as.integer(ploidy) if (ref_dosage) dose <- ploidy - dose #### Individuals in rows, markers in columns #### - base::rownames(dose) <- snp_names - geno <- base::t(dose) - if (base::is.null(base::rownames(geno)) || base::any(base::rownames(geno) == "")) + rownames(dose) <- snp_names + geno <- t(dose) + if (is.null(rownames(geno)) || any(rownames(geno) == "")) stop("Sample names could not be read from the VCF.") #### Call-rate filtering #### - n_samples_in <- base::nrow(geno) - n_snps_in <- base::ncol(geno) + n_samples_in <- nrow(geno) + n_snps_in <- ncol(geno) if (min_marker_call_rate > 0) { - marker_cr <- base::colMeans(!base::is.na(geno)) + marker_cr <- colMeans(!is.na(geno)) geno <- geno[, marker_cr >= min_marker_call_rate, drop = FALSE] - if (base::ncol(geno) == 0) + if (ncol(geno) == 0) stop("No markers passed min_marker_call_rate = ", min_marker_call_rate, ".") } if (min_sample_call_rate > 0) { - sample_cr <- base::rowMeans(!base::is.na(geno)) + sample_cr <- rowMeans(!is.na(geno)) geno <- geno[sample_cr >= min_sample_call_rate, , drop = FALSE] - if (base::nrow(geno) == 0) + if (nrow(geno) == 0) stop("No individuals passed min_sample_call_rate = ", min_sample_call_rate, ".") } #### Verbose report #### if (verbose) { - base::cat("\n=== VCF Conversion Report ===\n") - base::cat("Ploidy: ", ploidy, "\n") - base::cat("Markers read: ", n_markers_in, "\n") - base::cat("Markers returned: ", base::ncol(geno), "\n") - base::cat("Individuals returned: ", base::nrow(geno), " of ", n_samples_in, "\n", sep = "") - base::cat(base::sprintf("Missing calls: %.2f%%\n", - 100 * base::mean(base::is.na(geno)))) + cat("\n=== VCF Conversion Report ===\n") + cat(sprintf("Ploidy: %s\n", ploidy)) + cat(sprintf("Markers read: %s\n", n_markers_in)) + cat(sprintf("Markers returned: %d\n", ncol(geno))) + cat(sprintf("Individuals returned: %d of %s\n", nrow(geno), n_samples_in)) + cat(sprintf("Missing calls: %.2f%%\n", 100 * mean(is.na(geno)))) if (n_wrong_ploidy > 0) - base::cat("Calls set to missing for wrong ploidy:", n_wrong_ploidy, "\n") - base::cat("Coding: dosage of the ", - if (ref_dosage) "reference" else "alternate", " allele (0..", ploidy, ")\n", - sep = "") + cat(sprintf("Calls set to missing for wrong ploidy: %s\n", n_wrong_ploidy)) + cat(sprintf("Coding: dosage of the %s allele (0..%s)\n", + if (ref_dosage) "reference" else "alternate", ploidy)) } #### Return in the requested layout #### if (format == "matrix") { - base::attr(geno, "ploidy") <- ploidy + attr(geno, "ploidy") <- ploidy return(geno) } diff --git a/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/find_parentage.Rd b/man/find_parentage.Rd index b8b9b27..078bd13 100644 --- a/man/find_parentage.Rd +++ b/man/find_parentage.Rd @@ -84,9 +84,9 @@ for the other methods). Ties are then ranked like any other candidate, so \code{show_ties} has no effect. Progeny with fewer than \code{top_n} candidates have \code{NA} in the remaining ranks. Status refers to rank 1 only.} -\item{founders_file}{Character, optional. Path to a one-column file of -founder IDs. Only used when \code{progeny_file} is a pedigree: founders with -both parents unknown are skipped instead of searched. Default is NULL.} +\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: 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 7fd114a..1d0b2bc 100644 --- a/man/validate_pedigree.Rd +++ b/man/validate_pedigree.Rd @@ -28,8 +28,8 @@ with \code{\link[=vcf_to_dosage]{vcf_to_dosage()}} using \code{ploidy}; a path t data.table with an id column followed by marker columns coded as allele-B dosage (0, 1, ..., ploidy).} -\item{founders_file}{Character, optional. Path to a one-column file listing -founder IDs. Founders with both parents coded as 0 are left unchanged. +\item{founders_file}{Optional. Founder IDs as a path to a one-column file, a +character vector of ids, or a data.frame (first column). Founders with both parents coded as 0 are left unchanged. Defaults to NULL.} \item{trio_error_threshold}{Numeric. Maximum Mendelian error percentage to diff --git a/tests/.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-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) +})