diff --git a/DESCRIPTION b/DESCRIPTION index fbc9b5b..0ee9cae 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,76 +1,65 @@ Package: BIGr -Title: Breeding Insight Genomics Functions for Polyploid and Diploid Species -Version: 0.8.0 -Authors@R: c(person(given='Alexander M.', - family='Sandercock', - email='sandercock.alex@gmail.com', - role=c('cre','aut')), - person(given='Cristiane', - family='Taniguti', - role = 'aut'), - person(given='Josue', - family='Chinchilla-Vargas', - role='aut'), - person(given='Shufen', - family='Chen', - role='ctb'), - person(given='Manoj', - family='Sapkota', - role='ctb'), - person(given='Meng', - family='Lin', - role='ctb'), - person(given='Dongyan', - family='Zhao', - role='ctb'), - person('University', "of Florida", - role=c('cph'), - comment = "Breeding Insight")) +Title: Breeding Insight Genomics Functions for Polyploid and Diploid + Species +Version: 0.8.1 +Authors@R: c( + person("Alexander M.", "Sandercock", , "sandercock.alex@gmail.com", role = c("cre", "aut")), + person("Cristiane", "Taniguti", role = "aut"), + person("Josue", "Chinchilla-Vargas", role = "aut"), + person("Shufen", "Chen", role = "ctb"), + person("Manoj", "Sapkota", role = "ctb"), + person("Meng", "Lin", role = "ctb"), + person("Dongyan", "Zhao", role = "ctb"), + person("University", "of Florida", role = "cph", + comment = "Breeding Insight") + ) Maintainer: Alexander M. Sandercock Description: Functions developed within Breeding Insight to analyze diploid and polyploid breeding and genetic data. 'BIGr' provides the - ability to filter variant call format (VCF) files, extract single nucleotide polymorphisms (SNPs) - from diversity arrays technology missing allele discovery count (DArT MADC) files, - and manipulate genotype data for both diploid and polyploid species. It - also serves as the core dependency for the 'BIGapp' 'Shiny' app, which - provides a user-friendly interface for performing routine genotype - analysis tasks such as dosage calling, filtering, principal component analysis (PCA), - genome-wide association studies (GWAS), and - genomic prediction. For more details about the included 'breedTools' - functions, see Funkhouser et al. (2017) , and - the 'updog' output format, see Gerard et al. (2018) . + ability to filter variant call format (VCF) files, extract single + nucleotide polymorphisms (SNPs) from diversity arrays technology + missing allele discovery count (DArT MADC) files, and manipulate + genotype data for both diploid and polyploid species. It also serves + as the core dependency for the 'BIGapp' 'Shiny' app, which provides a + user-friendly interface for performing routine genotype analysis tasks + such as dosage calling, filtering, principal component analysis (PCA), + genome-wide association studies (GWAS), and genomic prediction. For + more details about the 'updog' output + format, see Gerard et al. (2018) . License: Apache License (>= 2) URL: https://github.com/Breeding-Insight/BIGr BugReports: https://github.com/Breeding-Insight/BIGr/issues -Encoding: UTF-8 -Roxygen: list(markdown = TRUE) -Depends: R (>= 4.4.0) -biocViews: +Depends: + R (>= 4.4.0) Imports: - parallel, + Biostrings, dplyr, + ggplot2, + parallel, + pwalign, Rdpack (>= 0.7), readr (>= 2.1.5), reshape2 (>= 1.4.4), rlang, - tidyr (>= 1.3.1), - vcfR (>= 1.15.0), Rsamtools, - Biostrings, - pwalign, - janitor, - quadprog, - tibble, + stats, stringr, - data.table + tibble, + tidyr (>= 1.3.1), + utils, + vcfR (>= 1.15.0) Suggests: + testthat (>= 3.0.0), covr, - ggplot2, spelling, - rmdformats, - knitr (>= 1.10), - rmarkdown, polyRAD, - testthat (>= 3.0.0) -RdMacros: Rdpack + rmdformats, + knitr, + rmarkdown +RdMacros: + Rdpack +biocViews: Config/roxygen2/version: 8.0.0 +Encoding: UTF-8 +Roxygen: list(markdown = TRUE) +RoxygenNote: 7.3.3 diff --git a/NAMESPACE b/NAMESPACE index 4d25dfe..4f200fe 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,17 +1,14 @@ # Generated by roxygen2: do not edit by hand -export(allele_freq_poly) export(calculate_Het) export(calculate_MAF) export(check_homozygous_trios) export(check_madc_sanity) -export(check_ped) export(check_replicates) export(dosage2vcf) export(dosage_ratios) export(filterMADC) export(filterVCF) -export(find_parentage) export(flip_dosage) export(get_countsMADC) export(imputation_concordance) @@ -20,13 +17,10 @@ export(madc2vcf_all) export(madc2vcf_multi) export(madc2vcf_targets) export(merge_MADCs) -export(solve_composition_poly) export(thinSNP) export(updog2vcf) -export(validate_pedigree) import(dplyr) import(parallel) -import(quadprog) import(stringr) import(tibble) import(tidyr) @@ -35,32 +29,16 @@ importFrom(Biostrings,DNAString) importFrom(Biostrings,reverseComplement) importFrom(Rdpack,reprompt) importFrom(Rsamtools,bgzip) -importFrom(data.table,CJ) -importFrom(data.table,as.data.table) -importFrom(data.table,copy) -importFrom(data.table,data.table) -importFrom(data.table,fread) -importFrom(data.table,is.data.table) -importFrom(data.table,rbindlist) -importFrom(data.table,set) importFrom(dplyr,"%>%") importFrom(dplyr,across) importFrom(dplyr,arrange) -importFrom(dplyr,bind_rows) importFrom(dplyr,case_when) -importFrom(dplyr,distinct) importFrom(dplyr,filter) -importFrom(dplyr,first) importFrom(dplyr,group_by) importFrom(dplyr,group_modify) -importFrom(dplyr,if_else) importFrom(dplyr,mutate) -importFrom(dplyr,n) -importFrom(dplyr,n_distinct) -importFrom(dplyr,row_number) importFrom(dplyr,select) importFrom(dplyr,summarise) -importFrom(dplyr,summarize) importFrom(dplyr,ungroup) importFrom(dplyr,where) importFrom(ggplot2,aes) @@ -70,14 +48,12 @@ importFrom(ggplot2,ggplot) importFrom(ggplot2,labs) importFrom(ggplot2,theme) importFrom(ggplot2,theme_minimal) -importFrom(janitor,clean_names) importFrom(pwalign,nucleotideSubstitutionMatrix) importFrom(pwalign,pairwiseAlignment) importFrom(readr,read_csv) importFrom(reshape2,dcast) importFrom(reshape2,melt) importFrom(rlang,sym) -importFrom(stats,cor) importFrom(stats,reorder) importFrom(stats,setNames) importFrom(tibble,as_tibble) diff --git a/NEWS.md b/NEWS.md index 309fc13..08e01b9 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,9 @@ +# BIGr 0.8.1 + +- Updated madc2vcf_all and madc2vcf_targets. Before, it was possible for POS to be exported as scientific notation instead of integers, and for negative POS values to be present for off target SNPs. POS are corrected to be integers, and SNPs with a negative POS value are removed. +- Remove BIGpopA functions - now it is a independent package: https://github.com/Breeding-Insight/BIGpopA +- Fixed `madc2vcf_all()` error "invalid substring arguments" that occurred with `add_others = TRUE` when an off-target ("Other") allele aligned to the reference with no mismatch positions remaining after the target SNP position was removed. The reference/alternate base lookups for off-target alleles are now guarded by the existing non-empty check, matching how the off-target Match alleles are already handled. + # BIGr 0.7.2 - Fixed manual text errors diff --git a/R/breedtools_functions.R b/R/breedtools_functions.R deleted file mode 100644 index 52a2d85..0000000 --- a/R/breedtools_functions.R +++ /dev/null @@ -1,182 +0,0 @@ -#' Compute allele frequencies for populations -#' -#' Computes allele frequencies for specified populations from SNP array data -#' coded as dosage of allele B. -#' -#' @param geno Numeric matrix of genotypes coded as dosage of allele B -#' \code{{0, 1, 2, ..., ploidy}}, with individuals in rows (named) and -#' SNPs in columns (named). -#' @param populations Named list of populations, each containing a character -#' vector of individual IDs belonging to that population. -#' @param ploidy Integer. Ploidy level of the species. Default is `2`. -#' -#' @return A data frame of allele frequencies with SNPs as rows and populations -#' as columns. -#' -#' @author Josué Chinchilla-Vargas -#' -#' @references Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP. -#' Estimation of genome-wide and locus-specific breed composition in pigs. -#' _Transl Anim Sci._ 2017;1(1):36–44. -#' -#' @examples -#' geno_matrix <- matrix( -#' c(4, 1, 4, 0, -#' 2, 2, 1, 3, -#' 0, 4, 0, 4, -#' 3, 3, 2, 2, -#' 1, 4, 2, 3), -#' nrow = 4, ncol = 5, byrow = FALSE, -#' dimnames = list(paste0("Ind", 1:4), paste0("S", 1:5)) -#' ) -#' pop_list <- list( -#' PopA = c("Ind1", "Ind2"), -#' PopB = c("Ind3", "Ind4") -#' ) -#' allele_freqs <- allele_freq_poly(geno = geno_matrix, populations = pop_list, ploidy = 4) -#' print(allele_freqs) -#' -#' @export -allele_freq_poly <- function(geno, populations, ploidy = 2) { - # Initialize result matrix - X <- matrix(NA, nrow = ncol(geno), ncol = length(populations)) - for (i in 1:length(populations)) { - pop_name <- names(populations[i]) - pop_geno <- geno[base::rownames(geno) %in% populations[[i]], ] - al_freq <- base::colMeans(pop_geno, na.rm = TRUE) / ploidy - X[, i] <- al_freq - } - base::colnames(X) <- base::names(populations) - base::rownames(X) <- base::colnames(geno) - return(X) -} - - -#' Solve breed composition for a single animal via quadratic programming -#' -#' Internal helper that solves the constrained OLS problem for one animal, -#' returning breed proportion estimates and the R² of the fit. -#' -#' @param Y Numeric vector of genotypes (named by SNP) for a single animal, -#' coded as dosage of allele B \code{{0, 1, 2, ..., ploidy}}. -#' @param X Numeric matrix of allele frequencies with SNPs as rows and -#' breeds/populations as columns. -#' -#' @return A named numeric vector of breed proportions plus `R2`. -#' -#' @references Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP. -#' Estimation of genome-wide and locus-specific breed composition in pigs. -#' _Transl Anim Sci._ 2017;1(1):36–44. -#' -#' @import quadprog -#' @importFrom stats cor -#' -#' @noRd -QPsolve <- function(Y, X) { - Ymod <- Y[!base::is.na(Y)] - Xmod <- X[base::names(Ymod), ] - p <- base::ncol(X) - nms <- base::colnames(X) - Rinv <- base::solve(base::chol(base::t(Xmod) %*% Xmod)) - C <- base::cbind(base::rep(1, p), base::diag(p)) - b2 <- c(1, base::rep(0, p)) - dd <- base::t(Ymod) %*% Xmod - qp <- quadprog::solve.QP(Dmat = Rinv, factorized = TRUE, dvec = dd, - Amat = C, bvec = b2, meq = 1) - beta <- qp$solution - rr <- stats::cor(Ymod, Xmod %*% beta)^2 - result <- c(beta, rr) - base::names(result) <- c(nms, "R2") - return(result) -} - - -#' Compute genome-wide breed composition -#' -#' Estimates genome-wide breed/ancestry composition for a batch of animals -#' using quadratic programming, with optional pedigree-assisted and -#' grouped-output modes. -#' -#' @param Y Numeric matrix of genotypes with individuals as rows and SNPs as -#' columns, coded as dosage of allele B \code{{0, 1, 2, ..., ploidy}}. -#' @param X Numeric matrix of allele frequencies with SNPs as rows and -#' breeds/populations as columns. -#' @param ped Optional data frame with pedigree information formatted with -#' columns `ID`, `Sire`, and `Dam`. When supplied, `QPsolve_par` is used -#' and only animals with genotyped parents are included. -#' @param groups Optional named list of IDs grouped by breed/population. -#' When supplied, results are returned as a named list partitioned by group. -#' @param mia Logical. Only applies when `ped` is supplied. If `TRUE`, returns -#' the inferred maternally inherited allele per locus per animal. Default `FALSE`. -#' @param sire Logical. Only applies when `ped` is supplied. If `TRUE`, returns -#' sire genotypes per locus per animal. Default `FALSE`. -#' @param dam Logical. Only applies when `ped` is supplied. If `TRUE`, returns -#' dam genotypes per locus per animal. Default `FALSE`. -#' @param ploidy Integer. Ploidy level of the species. Default is `2`. -#' -#' @return A data frame, or a named list of data frames when `groups` is -#' supplied, containing breed/ancestry composition estimates. -#' -#' @author Josué Chinchilla-Vargas -#' -#' @references Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP. -#' Estimation of genome-wide and locus-specific breed composition in pigs. -#' _Transl Anim Sci._ 2017;1(1):36–44. -#' -#' @examples -#' allele_freqs_matrix <- matrix( -#' c(0.625, 0.500, -#' 0.500, 0.500, -#' 0.500, 0.500, -#' 0.750, 0.500, -#' 0.625, 0.625), -#' nrow = 5, ncol = 2, byrow = TRUE, -#' dimnames = list(paste0("SNP", 1:5), c("VarA", "VarB")) -#' ) -#' val_geno_matrix <- matrix( -#' c(2, 1, 2, 3, 4, -#' 3, 4, 2, 3, 0), -#' nrow = 2, ncol = 5, byrow = TRUE, -#' dimnames = list(paste0("Test", 1:2), paste0("SNP", 1:5)) -#' ) -#' composition <- solve_composition_poly(Y = val_geno_matrix, -#' X = allele_freqs_matrix, -#' ploidy = 4) -#' print(composition) -#' -#' @import quadprog -#' -#' @export -solve_composition_poly <- function(Y, - X, - ped = NULL, - groups = NULL, - mia = FALSE, - sire = FALSE, - dam = FALSE, - ploidy = 2) { - # Transpose so Y is SNPs x animals, as required internally - Y <- base::t(Y) - # Retain only SNPs present in X - Y <- Y[base::rownames(Y) %in% base::rownames(X), ] - - if (!base::is.null(ped)) { - # Pedigree-assisted mode: use QPsolve_par for animals with genotyped parents - mat_results <- base::lapply(base::colnames(Y), QPsolve_par, Y, X, ped, - mia = mia, sire = sire, dam = dam) - mat_results_tab <- base::do.call(rbind, mat_results) - return(mat_results_tab) - - } else if (!base::is.null(groups)) { - # Grouped mode: standard computation, results partitioned by group - Y <- Y / ploidy - grouped_results <- base::lapply(groups, QPseparate, Y, X) - return(grouped_results) - - } else { - # Standard mode: unsegregated computation across all animals - Y <- Y / ploidy - results <- base::t(base::apply(Y, 2, QPsolve, X)) - return(results) - } -} diff --git a/R/check_ped.R b/R/check_ped.R deleted file mode 100644 index 259b696..0000000 --- a/R/check_ped.R +++ /dev/null @@ -1,293 +0,0 @@ -#' Check and Correct Common Pedigree Errors -#' -#' Reads a 3-column pedigree file (id, male_parent, female_parent) and performs -#' quality checks, optionally correcting detected errors. Exact duplicates and -#' missing parents are always corrected. Conflicting trios and inconsistent sex -#' roles are corrected when their respective arguments are TRUE. Cycles are -#' reported only and must be resolved manually. -#' -#' @param ped.file Path to the pedigree text file (TSV/CSV/TXT), OR a -#' data.frame / data.table with columns: id, male_parent, female_parent. -#' @param seed Optional integer seed for reproducibility. -#' @param verbose Logical. If TRUE (default), prints the report to the console. -#' @param correct_conflicting_trios Logical. If TRUE (default), sets conflicting -#' male_parent and female_parent to 0 and collapses to one row per ID. -#' @param correct_inconsistent_sex_roles Logical. If TRUE (default), sets -#' male_parent and female_parent to 0 for rows involving IDs found as both, -#' then removes any resulting exact duplicates. -#' -#' @return An invisible named list of data frames: -#' \describe{ -#' \item{exact_duplicates}{Exact duplicate rows found in the input.} -#' \item{conflicting_trios}{IDs with conflicting male_parent or female_parent assignments.} -#' \item{inconsistent_sex_roles}{Rows where a conflicting ID appears as male_parent or female_parent.} -#' \item{missing_parents}{Parent IDs absent from id, added as founders.} -#' \item{dependencies}{Cycles detected in the pedigree. Must be resolved manually.} -#' \item{corrected_pedigree}{Corrected pedigree table.} -#' } -#' -#' @examples -#' ped_file <- system.file("check_ped_test.txt", package = "BIGr") -#' ped_errors <- check_ped(ped.file = ped_file, seed = 101919, verbose = FALSE) -#' -#' # Also accepts a data.table directly -#' library(data.table) -#' ped_dt <- data.table(id = c("A","B","C"), -#' male_parent = c("0","0","A"), -#' female_parent = c("0","0","B")) -#' ped_errors <- check_ped(ped.file = ped_dt, verbose = FALSE) -#' -#' @author Josue Chinchilla-Vargas -#' -#' @importFrom dplyr %>% mutate filter group_by ungroup summarize distinct bind_rows select first n n_distinct if_else row_number -#' @importFrom stats setNames -#' @importFrom utils read.table -#' @importFrom janitor clean_names -#' @export -check_ped <- function(ped.file, - seed = NULL, - verbose = TRUE, - correct_conflicting_trios = TRUE, - correct_inconsistent_sex_roles = TRUE) { - - #### setup #### - if (!is.null(seed)) set.seed(seed) - - # Accept file path OR in-memory data.frame / data.table - if (is.character(ped.file) && length(ped.file) == 1 && file.exists(ped.file)) { - data <- utils::read.table(ped.file, header = TRUE) - data <- janitor::clean_names(data) - } else if (is.data.frame(ped.file) || data.table::is.data.table(ped.file)) { - data <- as.data.frame(ped.file) - data <- janitor::clean_names(data) - } else { - stop("ped.file must be a valid file path (character) or a data.frame / data.table.") - } - - required_cols <- c("id", "male_parent", "female_parent") - missing_cols <- setdiff(required_cols, colnames(data)) - if (length(missing_cols) > 0) { - stop( - "Input is missing required column(s): ", - paste(missing_cols, collapse = ", "), - ".\nExpected columns: id, male_parent, female_parent." - ) - } - extra_cols <- setdiff(names(data), required_cols) - data <- data[, c(required_cols, extra_cols)] - - data <- data %>% - dplyr::mutate( - id = as.character(id), - male_parent = as.character(male_parent), - female_parent = as.character(female_parent) - ) - - data <- data %>% dplyr::mutate(row_number = dplyr::row_number(), .before = id) - - errors <- list() - missing_parents <- data.frame( - row_number = integer(), - id = character(), - male_parent = character(), - female_parent = character(), - stringsAsFactors = FALSE - ) - - #### check 1: exact duplicates (always fixed) #### - exact_duplicates <- data[ - duplicated(data %>% dplyr::select(-row_number)) | - duplicated(data %>% dplyr::select(-row_number), fromLast = TRUE), - ] - if (nrow(exact_duplicates) > 0) { - data <- data %>% - dplyr::select(-row_number) %>% - dplyr::distinct() %>% - dplyr::mutate(row_number = dplyr::row_number(), .before = id) - } - - #### check 2: conflicting trios #### - repeated_ids <- data %>% - dplyr::group_by(id) %>% - dplyr::filter(dplyr::n() > 1) %>% - dplyr::ungroup() - - conflicting_trios_ids <- repeated_ids %>% - dplyr::group_by(id) %>% - dplyr::filter(dplyr::n_distinct(male_parent) > 1 | - dplyr::n_distinct(female_parent) > 1) %>% - dplyr::ungroup() - - if (correct_conflicting_trios && nrow(conflicting_trios_ids) > 0) { - data <- data %>% - dplyr::group_by(id) %>% - dplyr::summarize( - row_number = dplyr::first(row_number), - male_parent = if (dplyr::n_distinct(male_parent) > 1) "0" else dplyr::first(male_parent), - female_parent = if (dplyr::n_distinct(female_parent) > 1) "0" else dplyr::first(female_parent), - .groups = "drop" - ) %>% - dplyr::select(row_number, id, male_parent, female_parent) - } - - conflicting_trios <- conflicting_trios_ids - - #### check 3: missing parents (always fixed) #### - for (i in seq_len(nrow(data))) { - id <- data$id[i] - male_parent <- data$male_parent[i] - female_parent <- data$female_parent[i] - - if (male_parent != "0" && male_parent != id && !male_parent %in% data$id) { - missing_parents <- rbind( - missing_parents, - data.frame(row_number = data$row_number[i], id = male_parent, - male_parent = "0", female_parent = "0", - stringsAsFactors = FALSE) - ) - } - if (female_parent != "0" && female_parent != id && !female_parent %in% data$id) { - missing_parents <- rbind( - missing_parents, - data.frame(row_number = data$row_number[i], id = female_parent, - male_parent = "0", female_parent = "0", - stringsAsFactors = FALSE) - ) - } - - if (male_parent == id || female_parent == id) { - errors <- append(errors, paste("Dependency: Individual", id, - "cannot be its own parent")) - } - } - - missing_parents <- dplyr::distinct(missing_parents) - if (nrow(missing_parents) > 0) { - data <- dplyr::bind_rows(data, missing_parents) - } - - #### check 4: inconsistent sex roles #### - male_ids <- unique(data$male_parent[data$male_parent != "0"]) - female_ids <- unique(data$female_parent[data$female_parent != "0"]) - conflicting_sex_ids <- intersect(male_ids, female_ids) - - inconsistent_sex_roles <- data %>% - dplyr::filter(male_parent %in% conflicting_sex_ids | - female_parent %in% conflicting_sex_ids) - - if (correct_inconsistent_sex_roles && length(conflicting_sex_ids) > 0) { - data <- data %>% - dplyr::mutate( - male_parent = dplyr::if_else(male_parent %in% conflicting_sex_ids, "0", male_parent), - female_parent = dplyr::if_else(female_parent %in% conflicting_sex_ids, "0", female_parent) - ) %>% - dplyr::distinct(id, male_parent, female_parent, .keep_all = TRUE) - } - - #### check 5: dependencies (cycles) -- reported only #### - detect_all_cycles <- function(data) { - adj_list <- lapply(data$id, function(x) { - row <- data[data$id == x, ] - c(row$male_parent, row$female_parent) - }) - names(adj_list) <- data$id - - dfs <- function(node, visited, rec_stack, path) { - visited[node] <- TRUE - rec_stack[node] <- TRUE - path <- append(path, node) - cycles <- list() - - for (neighbor in adj_list[[node]]) { - if (neighbor %in% names(adj_list)) { - if (!visited[neighbor]) { - cycles <- append(cycles, dfs(neighbor, visited, rec_stack, path)) - } else if (rec_stack[neighbor]) { - cycle_start <- match(neighbor, path) - cycles <- append(cycles, list(path[cycle_start:length(path)])) - } - } - } - rec_stack[node] <- FALSE - return(cycles) - } - - visited <- stats::setNames(rep(FALSE, length(adj_list)), names(adj_list)) - rec_stack <- stats::setNames(rep(FALSE, length(adj_list)), names(adj_list)) - all_cycles <- list() - - for (node in names(adj_list)) { - if (!visited[node]) { - node_cycles <- dfs(node, visited, rec_stack, character()) - if (length(node_cycles) > 0) - all_cycles <- append(all_cycles, node_cycles) - } - } - return(all_cycles) - } - - cycles <- detect_all_cycles(data) - if (length(cycles) > 0) { - for (cycle_group in cycles) { - cycle_ids <- unique(unlist(cycle_group)) - errors <- append(errors, - paste("Cycle detected involving IDs:", - paste(cycle_ids, collapse = " -> "))) - } - } - - #### compile findings #### - input_ped_report <- list( - exact_duplicates = exact_duplicates, - conflicting_trios = conflicting_trios, - inconsistent_sex_roles = inconsistent_sex_roles, - missing_parents = missing_parents, - dependencies = data.frame(dependency = unique(unlist(errors)), - stringsAsFactors = FALSE), - corrected_pedigree = data %>% dplyr::select(-row_number) - ) - - #### output #### - if (verbose) { - cat("\n=== Pedigree Quality Check Report ===\n") - - if (nrow(exact_duplicates) > 0) { - cat("\nExact duplicate trios detected (removed in corrected pedigree):\n") - print(exact_duplicates) - } else cat("\nNo exact duplicate trios found.\n") - - if (nrow(conflicting_trios) > 0) { - cat("\nConflicting trios detected:\n") - print(conflicting_trios) - if (correct_conflicting_trios) { - cat(" -> parents set to 0 and collapsed to one row in corrected pedigree.\n") - } else { - cat(" -> correct_conflicting_trios = FALSE: left as-is in corrected pedigree.\n") - } - } else cat("\nNo conflicting trios found.\n") - - if (nrow(missing_parents) > 0) { - cat("\nParents missing as IDs (added as founders in corrected pedigree):\n") - print(missing_parents) - } else cat("\nNo missing parents found.\n") - - if (nrow(inconsistent_sex_roles) > 0) { - cat("\nIDs found as both male_parent and female_parent:\n") - print(inconsistent_sex_roles) - if (correct_inconsistent_sex_roles) { - cat(" -> parent fields set to 0 for conflicting IDs in corrected pedigree.\n") - } else { - cat(" -> correct_inconsistent_sex_roles = FALSE: left as-is.\n") - } - } else cat("\nNo IDs found as both male_parent and female_parent.\n") - - if (nrow(input_ped_report$dependencies) > 0) { - cat("\nDependencies detected (must be resolved manually):\n") - print(input_ped_report$dependencies) - } else cat("\nNo dependencies detected.\n") - - cat("\nThe corrected pedigree is included in the returned list as corrected_pedigree.\n") - } - - invisible(input_ped_report) -} diff --git a/R/find_parentage.R b/R/find_parentage.R deleted file mode 100644 index 758eb02..0000000 --- a/R/find_parentage.R +++ /dev/null @@ -1,443 +0,0 @@ -#' Find Parentage Assignments for Progeny -#' -#' Assigns the most likely parent(s) to each progeny from SNP genotype data -#' using Mendelian error rates or homozygous mismatch rates. Parents or progeny -#' absent from the genotype file are removed with a warning. -#' -#' @param genotypes_file Path to a TSV/CSV/TXT file, OR a data.frame / -#' data.table with an 'id' column followed by marker columns coded as 0, 1, 2. -#' @param parents_file Path to a TSV/CSV/TXT file, OR a data.frame / -#' data.table with an 'id' column and an optional 'sex' column -#' ('M', 'F', or 'A'). If absent, all parents are treated as ambiguous. -#' @param progeny_file Path to a TSV/CSV/TXT file, OR a data.frame / -#' data.table with an 'id' column. -#' @param method Character. One of \code{"best_male_parent"}, -#' \code{"best_female_parent"}, \code{"best_match"}, or -#' \code{"best_pair"} (default). -#' @param min_markers Integer. Minimum markers required; fewer flags -#' \code{low_markers} (default: \code{10}). -#' @param error_threshold Numeric. Maximum mismatch percentage; exceeded values -#' flag \code{high_error} (default: \code{5.0}). Must be between 0 and 100. -#' @param show_ties Logical. If \code{TRUE}, tied best pairs are appended as -#' suffix columns. Default is \code{TRUE}. -#' @param allow_parent_selfing Logical. If \code{FALSE}, candidate pairs with -#' identical male and female parent IDs are excluded. Applies only when -#' \code{method = "best_pair"}. Default is \code{FALSE}. -#' @param exclude_self_match Logical. If \code{TRUE}, each progeny ID is -#' excluded from its own candidate parent set, preventing self-matches when -#' progeny are also present in the parents file. Default is \code{TRUE}. -#' @param verbose Logical. If \code{TRUE}, prints progress and summary. -#' Default is \code{TRUE}. -#' @param plot_results Logical. If \code{TRUE}, plots the Mendelian error -#' distribution. Requires \code{ggplot2}. Default is \code{TRUE}. -#' -#' @return A named list (returned invisibly) with elements: -#' \describe{ -#' \item{pass}{Progeny with a confident parentage assignment.} -#' \item{high_error}{Progeny whose best assignment exceeds the error threshold.} -#' \item{low_markers}{Progeny with insufficient markers for a valid assignment.} -#' \item{full_results}{Complete data.table with all progeny and all output columns.} -#' \item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} -#' } -#' -#' @author Josue Chinchilla-Vargas -#' -#' @importFrom data.table fread copy CJ rbindlist set data.table as.data.table is.data.table -#' @export -find_parentage <- function(genotypes_file, parents_file, progeny_file, - method = "best_pair", - min_markers = 10, - error_threshold = 5.0, - show_ties = TRUE, - allow_parent_selfing = FALSE, - exclude_self_match = TRUE, - verbose = TRUE, - plot_results = TRUE) { - - ## silence R CMD check NOTEs - id <- sex <- male_parent <- female_parent <- NULL - mendelian_error_pct <- plot_status <- status <- NULL - - #### Input Validation #### - allowed_methods <- c("best_male_parent", "best_female_parent", - "best_match", "best_pair") - if (!method %in% allowed_methods) - stop("Method must be one of: ", paste(allowed_methods, collapse = ", ")) - if (min_markers < 1) - stop("min_markers must be a positive integer.") - if (error_threshold < 0 || error_threshold > 100) - stop("error_threshold must be between 0 and 100.") - - # Accept file path OR in-memory data.frame / data.table - read_flex <- function(x, label, ...) { - if (is.character(x) && length(x) == 1) { - if (!file.exists(x)) - stop("Error reading input files. Ensure paths are correct and files are TXT/TSV/CSV.") - data.table::fread(x, sep = "auto", ...) - } else if (is.data.frame(x) || data.table::is.data.table(x)) { - data.table::as.data.table(x) - } else { - stop(label, " must be a file path (character) or a data.frame / data.table.") - } - } - - tryCatch({ - genos <- read_flex(genotypes_file, "genotypes_file") - all_parents <- read_flex(parents_file, "parents_file") - progeny_candidates <- read_flex(progeny_file, "progeny_file") - }, error = function(e) { - stop("Error reading input files. Ensure paths are correct and files are TXT/CSV/TSV.") - }) - - valid_ids <- genos$id - removed_parents <- base::setdiff(all_parents$id, valid_ids) - if (base::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) { - 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)) { - 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], - by = "id") - female_parent_candidates <- base::unique(all_parents[sex %in% c("F", "A", "NA")], - by = "id") - - if (base::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 && - 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) - stop("No valid progeny candidates remain after filtering.") - - #### Pre-compute genotype matrices once #### - genos_mat <- base::as.matrix(genos, rownames = "id") - genos_hom <- data.table::copy(genos) - marker_cols <- base::setdiff(base::names(genos_hom), "id") - for (col in marker_cols) - genos_hom[base::get(col) == 1, (col) := NA_integer_] - genos_hom_mat <- base::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")) - } - - #### Logic for Homozygous Matching Methods #### - if (method %in% c("best_male_parent", "best_female_parent", "best_match")) { - parent_ids <- base::switch(method, - "best_male_parent" = male_parent_candidates$id, - "best_female_parent" = female_parent_candidates$id, - "best_match" = base::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] - - n_progeny <- base::nrow(progeny_genos) - results_dt <- data.table::data.table( - id = base::rownames(progeny_genos), - best_match = NA_character_, - mendelian_error_pct = NA_real_, - markers_tested = NA_integer_, - status = NA_character_ - ) - - for (i in base::seq_len(n_progeny)) { - prog_id <- base::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), - 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)) - percent_mismatch <- (mismatches / comparisons) * 100 - percent_mismatch[base::is.nan(percent_mismatch)] <- NA - - if (exclude_self_match) { - self_idx <- base::rownames(parent_genos) == prog_id - percent_mismatch[self_idx] <- NA_real_ - } - - best_idx <- base::which.min(percent_mismatch) - if (base::length(best_idx) == 0) { - data.table::set(results_dt, i, "markers_tested", 0L) - data.table::set(results_dt, i, "status", "low_markers") - 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]) - data.table::set(results_dt, i, "mendelian_error_pct", best_error) - data.table::set(results_dt, i, "markers_tested", base::as.integer(best_markers)) - data.table::set(results_dt, i, "status", assign_status(best_markers, best_error)) - } - final_df <- results_dt - } - - #### Logic for Best Pair Method #### - if (method == "best_pair") { - parent_pairs <- data.table::CJ(male_parent = male_parent_candidates$id, - 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 (base::nrow(parent_pairs) == 0) stop("No valid parent pairs to test.") - - male_parent_genos_mat <- genos_mat[parent_pairs$male_parent, , drop = FALSE] - female_parent_genos_mat <- genos_mat[parent_pairs$female_parent, , drop = FALSE] - progeny_ids <- progeny_candidates$id - progeny_mat <- genos_mat[progeny_ids, , drop = FALSE] - n_progeny <- base::nrow(progeny_mat) - n_pairs <- base::nrow(parent_pairs) - - mismatch_mat <- base::matrix( - base::vapply(base::seq_len(n_progeny), function(j) { - progeny_vec <- progeny_mat[j, ] - progeny_pair_mat <- base::matrix(progeny_vec, - nrow = n_pairs, - ncol = base::ncol(progeny_mat), - byrow = TRUE) - base::rowSums( - (male_parent_genos_mat == 0 & female_parent_genos_mat == 0 & progeny_pair_mat > 0) | - (male_parent_genos_mat == 2 & female_parent_genos_mat == 2 & progeny_pair_mat < 2) | - ((male_parent_genos_mat == 0 & female_parent_genos_mat == 1) | - (male_parent_genos_mat == 1 & female_parent_genos_mat == 0)) & (progeny_pair_mat == 2) | - ((male_parent_genos_mat == 2 & female_parent_genos_mat == 1) | - (male_parent_genos_mat == 1 & female_parent_genos_mat == 2)) & (progeny_pair_mat == 0) | - ((male_parent_genos_mat == 0 & female_parent_genos_mat == 2) | - (male_parent_genos_mat == 2 & female_parent_genos_mat == 0)) & (progeny_pair_mat != 1), - na.rm = TRUE - ) - }, numeric(n_pairs)), - nrow = n_pairs, ncol = n_progeny - ) - - comparison_mat <- base::matrix( - base::vapply(base::seq_len(n_progeny), function(j) { - progeny_vec <- progeny_mat[j, ] - progeny_pair_mat <- base::matrix(progeny_vec, - nrow = n_pairs, - ncol = base::ncol(progeny_mat), - byrow = TRUE) - base::rowSums(!base::is.na(male_parent_genos_mat) & - !base::is.na(female_parent_genos_mat) & - !base::is.na(progeny_pair_mat)) - }, numeric(n_pairs)), - nrow = n_pairs, ncol = n_progeny - ) - - pct_mismatch_mat <- (mismatch_mat / comparison_mat) * 100 - pct_mismatch_mat[base::is.nan(pct_mismatch_mat)] <- NA - - results_dt <- data.table::data.table( - id = progeny_ids, - male_parent = NA_character_, - female_parent = NA_character_, - mendelian_error_pct = NA_character_, - markers_tested = NA_integer_, - status = NA_character_ - ) - - results_list <- base::vector("list", n_progeny) - for (j in base::seq_len(n_progeny)) { - prog_id <- progeny_ids[j] - percent_mismatch <- pct_mismatch_mat[, j] - comparisons <- comparison_mat[, j] - - if (exclude_self_match) { - self_pair_idx <- parent_pairs$male_parent == prog_id | - parent_pairs$female_parent == prog_id - percent_mismatch[self_pair_idx] <- NA_real_ - } - - min_mismatch_val <- base::min(percent_mismatch, na.rm = TRUE) - - if (base::is.infinite(min_mismatch_val)) { - data.table::set(results_dt, j, "markers_tested", 0L) - data.table::set(results_dt, j, "status", "low_markers") - next - } - - best_indices <- base::which(percent_mismatch == min_mismatch_val) - if (base::length(best_indices) > 1) { - best_markers_per_pair <- comparisons[best_indices] - max_markers <- base::max(best_markers_per_pair) - best_indices <- best_indices[best_markers_per_pair == max_markers] - } - - best_pairs <- parent_pairs[best_indices] - best_markers <- comparisons[best_indices[1]] - best_error <- base::round(min_mismatch_val, 2) - a_status <- assign_status(best_markers, best_error) - - if (!show_ties && base::nrow(best_pairs) > 1) { - warning("Progeny '", prog_id, "' has ", base::nrow(best_pairs), - " tied best pairs. Only one is reported as show_ties=FALSE.", - call. = FALSE) - } - - num_to_report <- base::min(base::nrow(best_pairs), - if (show_ties) base::nrow(best_pairs) else 1) - - data.table::set(results_dt, j, "male_parent", best_pairs$male_parent[1]) - data.table::set(results_dt, j, "female_parent", best_pairs$female_parent[1]) - data.table::set(results_dt, j, "mendelian_error_pct", base::sprintf("%.2f", min_mismatch_val)) - data.table::set(results_dt, j, "markers_tested", base::as.integer(best_markers)) - data.table::set(results_dt, j, "status", a_status) - - if (show_ties && num_to_report > 1) { - tie_row <- base::list(id = prog_id) - for (k in base::seq(2, num_to_report)) { - tie_row[[base::paste0("male_parent_", k)]] <- best_pairs$male_parent[k] - tie_row[[base::paste0("female_parent_", k)]] <- best_pairs$female_parent[k] - tie_row[[base::paste0("mendelian_error_pct_", k)]] <- min_mismatch_val - tie_row[[base::paste0("markers_tested_", k)]] <- comparisons[best_indices[k]] - } - results_list[[j]] <- data.table::as.data.table(tie_row) - } - } - - tie_rows <- data.table::rbindlist( - base::Filter(Negate(base::is.null), results_list), - fill = TRUE, - use.names = TRUE - ) - if (base::nrow(tie_rows) > 0) { - final_df <- merge(results_dt, tie_rows, by = "id", all.x = TRUE) - for (col in base::names(final_df)) - data.table::set(final_df, which(final_df[[col]] == ""), col, NA_character_) - } else { - final_df <- results_dt - } - } - - #### Compile named list #### - output_list <- list( - pass = final_df[status == "pass"], - high_error = final_df[status == "high_error"], - low_markers = final_df[status == "low_markers"], - full_results = final_df, - plot = NULL - ) - - #### Verbose output #### - if (verbose) { - total_progeny <- base::nrow(final_df) - base::cat("\n=== Parentage Assignment Report ===\n") - base::cat("\nTotal progeny evaluated:", total_progeny, "\n") - base::cat("Method:", method, " | ", - "Error threshold:", error_threshold, "% | ", - "Minimum markers:", min_markers, "\n") - - n_pass <- base::nrow(output_list$pass) - if (n_pass > 0) { - base::cat(base::sprintf("\n%d progeny passed (%.1f%%).\n", - n_pass, (n_pass / total_progeny) * 100)) - } else { - base::cat("\nNo progeny passed.\n") - } - - n_high <- base::nrow(output_list$high_error) - if (n_high > 0) { - base::cat(base::sprintf("\n%d progeny flagged high_error (%.1f%%):\n", - n_high, (n_high / total_progeny) * 100)) - base::print(output_list$high_error) - } else { - base::cat("\nNo progeny flagged for high error.\n") - } - - n_low <- base::nrow(output_list$low_markers) - if (n_low > 0) { - base::cat(base::sprintf("\n%d progeny flagged low_markers (%.1f%%):\n", - n_low, (n_low / total_progeny) * 100)) - base::print(output_list$low_markers) - } else { - base::cat("\nNo progeny flagged for low marker count.\n") - } - - base::cat("\nFull results are included in the returned list as $full_results.\n") - } - - #### Plot Results #### - if (plot_results) { - 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$mendelian_error_pct)] - plot_df$mendelian_error_pct <- base::as.numeric(plot_df$mendelian_error_pct) - plot_df$plot_status <- base::ifelse( - plot_df$status == "pass", "pass", - base::ifelse( - plot_df$status == "high_error", "high_error", - base::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) - - threshold_label <- base::paste0( - "Error Threshold: ", error_threshold, "% | ", - "Pass: ", n_pass, " | ", - "High Error: ", n_high, " | ", - "Low Markers: ", n_low - ) - - p <- ggplot2::ggplot( - plot_df, - ggplot2::aes(x = mendelian_error_pct, fill = plot_status) - ) + - ggplot2::geom_histogram(binwidth = 1, color = "white", alpha = 0.9) + - ggplot2::geom_vline(xintercept = error_threshold, - linetype = "dashed", color = "black", linewidth = 1) + - ggplot2::scale_x_continuous(breaks = seq(0, 100, by = 5)) + - ggplot2::scale_y_continuous(breaks = seq(0, 10000, by = 5)) + - ggplot2::scale_fill_manual( - values = c("pass" = "#339900", - "high_error" = "#cc3333", - "low_markers" = "#F1C40F", - "other" = "#BDC3C7"), - labels = c("pass" = "Pass", - "high_error" = "High Error", - "low_markers" = "Low Markers", - "other" = "Other") - ) + - ggplot2::labs( - title = "Parentage Mendelian Error Distribution", - subtitle = base::paste0("Progeny Tested: ", n_total, - "\n \n", threshold_label), - x = "Mendelian Error (%)", - y = "Number of Progeny", - fill = "Status" - ) + - ggplot2::theme_classic(base_size = 13) + - ggplot2::theme(legend.position = "top") - - base::print(p) - output_list$plot <- p - } - } - - return(base::invisible(output_list)) -} diff --git a/R/madc2vcf_all.R b/R/madc2vcf_all.R index e223977..b11ac9e 100644 --- a/R/madc2vcf_all.R +++ b/R/madc2vcf_all.R @@ -652,10 +652,11 @@ compare <- function(one_tag, botloci, alignment_score_thr = 40, mi_df = NULL, ad pos_alt_idx <- pos_alt_idx[-rm_target_other] } } - other_ref_base <- substring(ref_seq, pos_ref_idx, pos_ref_idx) - other_alt_base <- substring(others_seq[j,]$AlleleSequence, pos_alt_idx, pos_alt_idx) # Cases found where the AltMatch is another alternative for the target SNP - they are discarted if(length(pos_ref_idx) >0){ + # Compute bases only when mismatch positions remain; substring() errors on integer(0) indices + other_ref_base <- substring(ref_seq, pos_ref_idx, pos_ref_idx) + other_alt_base <- substring(others_seq[j,]$AlleleSequence, pos_alt_idx, pos_alt_idx) # If Match sequences have N, do not consider as polymorphism if(any(!other_alt_base %in% c("A", "T", "C", "G"))) { other_ref_base <- other_ref_base[-which(!other_alt_base %in% c("A", "T", "C", "G"))] @@ -855,7 +856,20 @@ create_VCF_body <- function(csv, vcf_body_new <- vcf_body_new[,-1] colnames(vcf_body_new) <- c("#CHROM", "POS", "ID", "REF", "ALT","QUAL", "FILTER", "INFO","FORMAT", colnames(csv)[-c(1:7)]) - vcf_body_new <- vcf_body_new[order(vcf_body_new[,1], vcf_body_new[,2]),] + + #Correct scientific notation POS and remove rows with negative POS information + pos_num <- suppressWarnings(as.numeric(vcf_body_new$POS)) + bad_pos <- is.na(pos_num) | !is.finite(pos_num) | pos_num < 1 | pos_num != floor(pos_num) + + if (any(bad_pos)) { + vmsg("Removed %s VCF records with invalid POS values.", verbose = verbose, level = 1, type = ">>", sum(bad_pos)) + vcf_body_new <- vcf_body_new[!bad_pos, , drop = FALSE] + pos_num <- pos_num[!bad_pos] + } + + ord <- order(vcf_body_new[[1]], pos_num) + vcf_body_new <- vcf_body_new[ord, , drop = FALSE] + vcf_body_new$POS <- sprintf("%.0f", pos_num[ord]) return(vcf_body_new) } diff --git a/R/madc2vcf_targets.R b/R/madc2vcf_targets.R index 81e0643..c126d08 100644 --- a/R/madc2vcf_targets.R +++ b/R/madc2vcf_targets.R @@ -128,6 +128,7 @@ #' @import dplyr #' @import tidyr #' @import tibble +#' @importFrom stats setNames #' @importFrom reshape2 melt dcast #' @importFrom utils write.table #' @importFrom Biostrings DNAString reverseComplement @@ -521,8 +522,19 @@ madc2vcf_targets <- function(madc_file, # Add # to the CHROM column name colnames(vcf_df)[1] <- "#CHROM" - # Sort - vcf_df <- vcf_df[order(vcf_df[,1],as.numeric(as.character(vcf_df[,2]))),] + # Sort and correct any potential scientific notation + pos_num <- suppressWarnings(as.numeric(vcf_df$POS)) + bad_pos <- is.na(pos_num) | !is.finite(pos_num) | pos_num < 1 | pos_num != floor(pos_num) + + if (any(bad_pos)) { + vmsg("Removed %s VCF records with invalid POS values.", verbose = verbose, level = 1, type = ">>", sum(bad_pos)) + vcf_df <- vcf_df[!bad_pos, , drop = FALSE] + pos_num <- pos_num[!bad_pos] + } + + ord <- order(vcf_df[[1]], pos_num) + vcf_df <- vcf_df[ord, , drop = FALSE] + vcf_df$POS <- sprintf("%.0f", pos_num[ord]) # Remove markers with NA CHROM/POS (unmatched in markers_info, Case 3) na_coord <- is.na(vcf_df[, 1]) | is.na(vcf_df$POS) diff --git a/R/validate_pedigree.R b/R/validate_pedigree.R deleted file mode 100644 index fc4cd0f..0000000 --- a/R/validate_pedigree.R +++ /dev/null @@ -1,431 +0,0 @@ -#' Validate Pedigree Trios Using Mendelian Error Analysis -#' -#' Validates parent-offspring trios against SNP genotype data using Mendelian -#' error rates. Identifies incorrect parentage assignments, suggests -#' best-matching replacements, and outputs a corrected pedigree. Founder trios -#' (both parents coded as 0) are preserved unchanged if a founders file is -#' supplied. Trios absent from the genotype file are retained as -#' no_genotype_data. -#' -#' @param pedigree_file Path to the pedigree file (TSV/CSV/TXT), OR a -#' data.frame / data.table with columns: id, male_parent, female_parent. -#' @param genotypes_file Path to the genotypes file (TSV/CSV/TXT), OR a -#' data.frame / data.table with an id column followed by marker columns -#' coded as 0, 1, 2. -#' @param founders_file Character, optional. Path to a one-column file listing -#' founder IDs. 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. -#' @param min_markers Integer. Minimum non-missing markers required to evaluate -#' a trio (default: 10). -#' @param single_parent_error_threshold Numeric. Maximum homozygous-marker -#' mismatch percentage for a parent to be considered acceptable -#' (default: 2.0). Must be between 0 and 100. -#' @param verbose Logical. If TRUE, prints progress, summary, and results to -#' the console (default: TRUE). -#' @param plot_results Logical. If TRUE, prints a histogram of trio Mendelian -#' error percentages with a threshold line (default: TRUE). -#' -#' @return An invisible named list with the following elements: -#' \describe{ -#' \item{pass}{Trios that passed the Mendelian error threshold.} -#' \item{fail}{Trios that failed the Mendelian error threshold.} -#' \item{low_markers}{Trios with insufficient markers for evaluation.} -#' \item{no_genotype_data}{Trios absent from the genotype file.} -#' \item{founders}{Trios identified as founders.} -#' \item{missing_parents}{Trios with one or both parents coded as 0 (non-founders).} -#' \item{full_results}{Complete data.table with all trios and all output columns.} -#' \item{corrected_pedigree}{Pedigree table after applying recommended corrections.} -#' \item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} -#' } -#' -#' @author Josue Chinchilla-Vargas -#' -#' @importFrom data.table fread copy data.table set rbindlist as.data.table is.data.table -#' @export -validate_pedigree <- function(pedigree_file, genotypes_file, - founders_file = NULL, - trio_error_threshold = 5.0, - min_markers = 10, - single_parent_error_threshold = 2.0, - verbose = TRUE, - plot_results = TRUE) { - - ## silence R CMD check NOTEs - id <- male_parent <- female_parent <- status <- trio_mendelian_error_pct <- NULL - plot_status <- recommended_correction <- NULL - - #### Input validation #### - if (trio_error_threshold < 0 || trio_error_threshold > 100) - stop("trio_error_threshold must be between 0 and 100") - if (single_parent_error_threshold < 0 || single_parent_error_threshold > 100) - stop("single_parent_error_threshold must be between 0 and 100") - - # Accept file path OR in-memory data.frame / data.table - read_flex <- function(x, label, ...) { - if (is.character(x) && length(x) == 1) { - if (!file.exists(x)) - stop("Error reading input files. Ensure paths are correct and files are TXT/TSV/CSV.") - data.table::fread(x, sep = "auto", ...) - } else if (is.data.frame(x) || data.table::is.data.table(x)) { - data.table::as.data.table(x) - } else { - stop(label, " must be a file path (character) or a data.frame / data.table.") - } - } - - tryCatch({ - pedigree <- read_flex(pedigree_file, "pedigree_file", colClasses = "character") - genos <- read_flex(genotypes_file, "genotypes_file") - }, error = function(e) { - stop("Error reading input files. Ensure paths are correct and files are TXT/TSV/CSV.") - }) - - #### Check required columns #### - 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) - stop("Pedigree file missing required columns: ", - base::paste(missing_cols, collapse = ", ")) - if (!"id" %in% base::names(genos)) - stop("Genotypes file must have an 'id' column") - - pedigree[, male_parent := as.character(male_parent)] - pedigree[, female_parent := as.character(female_parent)] - 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) - } - - #### Build genotype matrices #### - genos_mat <- base::as.matrix(genos, rownames = "id") - genos_hom <- data.table::copy(genos) - marker_cols <- base::setdiff(base::names(genos_hom), "id") - for (col in marker_cols) - genos_hom[base::get(col) == 1, (col) := NA_integer_] - genos_hom_mat <- base::as.matrix(genos_hom, rownames = "id") - - #### Identify trios missing from the genotype file #### - valid_ids <- as.character(genos$id) - has_geno <- pedigree[id %in% valid_ids & - (male_parent %in% valid_ids | male_parent == "0") & - (female_parent %in% valid_ids | female_parent == "0")] - 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") - pedigree <- has_geno - if (base::nrow(pedigree) == 0) - stop("No valid trios remain after filtering for genotype availability.") - - #### Find best matching parent via homozygous mismatch #### - find_best_parent <- function(prog_id, exclude_ids = base::character(0)) { - candidates <- base::setdiff(base::rownames(genos_hom_mat), - c(prog_id, exclude_ids)) - if (base::length(candidates) == 0) - return(base::list(id = NA_character_, error_pct = NA_real_)) - prog_hom <- genos_hom_mat[prog_id, ] - errors <- base::sapply(candidates, function(cand_id) { - cand_hom <- genos_hom_mat[cand_id, ] - comparisons <- base::sum(!base::is.na(cand_hom) & !base::is.na(prog_hom)) - if (comparisons == 0) return(NA_real_) - (base::sum(cand_hom != prog_hom, na.rm = TRUE) / comparisons) * 100 - }) - if (base::all(base::is.na(errors))) - return(base::list(id = NA_character_, error_pct = NA_real_)) - best_idx <- base::which.min(errors) - base::list(id = candidates[best_idx], - error_pct = base::round(errors[best_idx], 2)) - } - - #### Main trio evaluation loop #### - results_list <- base::lapply(base::seq_len(base::nrow(pedigree)), function(i) { - prog_id <- pedigree$id[i] - male_parent_id <- pedigree$male_parent[i] - female_parent_id <- pedigree$female_parent[i] - correction_decision <- "none" - error_pct <- NA_real_ - status <- "no_data" - markers_tested <- 0L - male_parent_error_pct <- NA_real_ - female_parent_error_pct <- NA_real_ - best_male_parent <- NA_character_ - best_male_parent_pct <- NA_real_ - best_female_parent <- NA_character_ - best_female_parent_pct <- NA_real_ - - if (male_parent_id == "0" && female_parent_id == "0" && - prog_id %in% founder_ids) { - status <- "founders" - correction_decision <- "none" - } else { - if (male_parent_id == "0" && female_parent_id == "0") { - status <- "missing_both_parents" - correction_decision <- "none" - best_m <- find_best_parent(prog_id, exclude_ids = character(0)) - best_male_parent <- best_m$id - best_male_parent_pct <- best_m$error_pct - best_f <- find_best_parent(prog_id, exclude_ids = c(best_m$id)) - best_female_parent <- best_f$id - best_female_parent_pct <- best_f$error_pct - } else if (male_parent_id == "0" && female_parent_id != "0") { - status <- "missing_male_parent" - correction_decision <- "none" - best_m <- find_best_parent(prog_id, exclude_ids = c(female_parent_id)) - best_male_parent <- best_m$id - best_male_parent_pct <- best_m$error_pct - } else if (male_parent_id != "0" && female_parent_id == "0") { - status <- "missing_female_parent" - correction_decision <- "none" - best_f <- find_best_parent(prog_id, exclude_ids = c(male_parent_id)) - best_female_parent <- best_f$id - best_female_parent_pct <- best_f$error_pct - } else { - progeny_vec <- genos_mat[prog_id, ] - male_parent_vec <- genos_mat[male_parent_id, ] - female_parent_vec <- genos_mat[female_parent_id, ] - mismatches <- base::sum( - (male_parent_vec == 0 & female_parent_vec == 0 & progeny_vec > 0) | - (male_parent_vec == 2 & female_parent_vec == 2 & progeny_vec < 2) | - ((male_parent_vec == 0 & female_parent_vec == 1) | - (male_parent_vec == 1 & female_parent_vec == 0)) & (progeny_vec == 2) | - ((male_parent_vec == 2 & female_parent_vec == 1) | - (male_parent_vec == 1 & female_parent_vec == 2)) & (progeny_vec == 0) | - ((male_parent_vec == 0 & female_parent_vec == 2) | - (male_parent_vec == 2 & female_parent_vec == 0)) & (progeny_vec != 1), - na.rm = TRUE - ) - markers_tested <- base::sum(!base::is.na(male_parent_vec) & - !base::is.na(female_parent_vec) & - !base::is.na(progeny_vec)) - if (markers_tested == 0) { - status <- "no_data" - correction_decision <- "none" - } else { - error_pct <- (mismatches / markers_tested) * 100 - if (markers_tested < min_markers) { - status <- "low_markers" - } else if (error_pct <= trio_error_threshold) { - status <- "pass" - correction_decision <- "none" - } else { - status <- "fail" - } - if (status %in% c("fail", "low_markers")) { - 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_parent_error_pct <- if (male_comparisons == 0) NA_real_ else - base::round((base::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_parent_error_pct <- if (female_comparisons == 0) NA_real_ else - base::round((base::sum(female_parent_hom != progeny_hom, na.rm = TRUE) / - female_comparisons) * 100, 2) - male_acceptable <- !is.na(male_parent_error_pct) && - male_parent_error_pct <= single_parent_error_threshold - female_acceptable <- !is.na(female_parent_error_pct) && - female_parent_error_pct <= single_parent_error_threshold - if (male_acceptable && female_acceptable) { - correction_decision <- "keep_both" - } else if (male_acceptable && !female_acceptable) { - correction_decision <- "remove_female_parent" - best_f <- find_best_parent(prog_id, exclude_ids = c(male_parent_id)) - best_female_parent <- best_f$id - best_female_parent_pct <- best_f$error_pct - } else if (!male_acceptable && female_acceptable) { - correction_decision <- "remove_male_parent" - best_m <- find_best_parent(prog_id, exclude_ids = c(female_parent_id)) - best_male_parent <- best_m$id - best_male_parent_pct <- best_m$error_pct - } else { - correction_decision <- "remove_both" - best_m <- find_best_parent(prog_id, exclude_ids = character(0)) - best_male_parent <- best_m$id - best_male_parent_pct <- best_m$error_pct - best_f <- find_best_parent(prog_id, exclude_ids = c(best_m$id)) - best_female_parent <- best_f$id - best_female_parent_pct <- best_f$error_pct - } - if (status == "low_markers") - correction_decision <- paste0("low_markers_", correction_decision) - } - } - } - } - - data.table::data.table( - 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_markers_tested = markers_tested, - status = status, - recommended_correction = correction_decision, - male_parent_hom_error_pct = male_parent_error_pct, - female_parent_hom_error_pct = female_parent_error_pct, - best_male_candidate = best_male_parent, - best_male_candidate_error_pct = best_male_parent_pct, - best_female_candidate = best_female_parent, - best_female_candidate_error_pct = best_female_parent_pct - ) - }) - - final_df <- data.table::rbindlist(results_list) - - #### Append no_genotype_data rows #### - if (base::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, - orig_female_parent = no_geno_rows$female_parent, - trio_mendelian_error_pct = NA_real_, - trio_markers_tested = 0L, - status = "no_genotype_data", - recommended_correction = "none", - male_parent_hom_error_pct = NA_real_, - female_parent_hom_error_pct = NA_real_, - best_male_candidate = NA_character_, - best_male_candidate_error_pct = NA_real_, - best_female_candidate = NA_character_, - best_female_candidate_error_pct = NA_real_ - ) - final_df <- data.table::rbindlist(list(final_df, no_geno_df)) - } - - #### Build corrected pedigree #### - corrected_pedigree <- data.table::copy(original_pedigree) - for (i in base::seq_len(base::nrow(final_df))) { - prog_id <- final_df$id[i] - decision <- final_df$recommended_correction[i] - row_idx <- base::which(corrected_pedigree$id == prog_id) - if (decision == "remove_male_parent") { - data.table::set(corrected_pedigree, row_idx, "male_parent", "0") - } else if (decision == "remove_female_parent") { - data.table::set(corrected_pedigree, row_idx, "female_parent", "0") - } else if (decision %in% c("remove_both", - "low_markers_remove_both", - "low_markers_remove_male_parent", - "low_markers_remove_female_parent")) { - if (grepl("male", decision)) - data.table::set(corrected_pedigree, row_idx, "male_parent", "0") - if (grepl("female", decision)) - data.table::set(corrected_pedigree, row_idx, "female_parent", "0") - if (decision %in% c("low_markers_remove_both", "remove_both")) { - data.table::set(corrected_pedigree, row_idx, "male_parent", "0") - data.table::set(corrected_pedigree, row_idx, "female_parent", "0") - } - } - } - - #### Summary output #### - if (verbose) { - total_trios <- base::nrow(final_df) - status_counts <- base::table(final_df$status) - base::cat("\n--- Trio Validation Summary ---\n") - base::cat("Total trios in pedigree:", total_trios, "\n") - for (s in base::names(status_counts)) - base::cat(base::sprintf("%-28s: %d (%.1f%%)\n", s, - status_counts[s], - (status_counts[s] / total_trios) * 100)) - base::cat("Error threshold:", trio_error_threshold, "%\n") - base::cat("Homozygous threshold:", single_parent_error_threshold, "%\n") - base::cat("Minimum markers required:", min_markers, "\n\n") - corrections <- base::table(final_df$recommended_correction) - base::cat("Correction summary:\n") - for (decision in base::names(corrections)) - if (decision != "none") - base::cat(" ", decision, ":", corrections[decision], "\n") - base::cat("\n") - base::print(final_df) - } - - #### Plot results #### - p <- NULL - if (plot_results) { - if (!requireNamespace("ggplot2", quietly = TRUE)) { - warning("ggplot2 is required for plot_results = TRUE. Please install it.", call. = FALSE) - } else { - plot_df <- final_df[!is.na(final_df$trio_mendelian_error_pct)] - plot_df$plot_status <- dplyr::case_when( - plot_df$recommended_correction %in% c("none", "keep_both", - "low_markers_keep_both") ~ "pass", - plot_df$recommended_correction %in% c("remove_male_parent", - "remove_female_parent", - "low_markers_remove_male_parent", - "low_markers_remove_female_parent") ~ "fail_one_parent", - plot_df$recommended_correction %in% c("remove_both", - "low_markers_remove_both") ~ "fail_both_parents", - TRUE ~ "other" - ) - 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) - threshold_label <- paste0( - "Mendelian Error Threshold: ", trio_error_threshold, "% | ", - "Lost: ", n_fail, " trios | ", - "Kept: ", n_pass, " trios" - ) - p <- ggplot2::ggplot(plot_df, - ggplot2::aes(x = trio_mendelian_error_pct, - fill = plot_status)) + - ggplot2::geom_histogram(binwidth = 1, color = "white", alpha = 0.9) + - ggplot2::geom_vline(xintercept = trio_error_threshold, - linetype = "dashed", color = "black", linewidth = 1) + - ggplot2::scale_x_continuous(breaks = seq(0, 100, by = 5)) + - ggplot2::scale_y_continuous(breaks = seq(0, 100, by = 5)) + - ggplot2::scale_fill_manual( - values = c("pass" = "#339900", - "fail_one_parent" = "#F1C40F", - "fail_both_parents" = "#cc3333", - "other" = "#BDC3C7"), - labels = c("pass" = "Pass", - "fail_one_parent" = "Fail - One Parent", - "fail_both_parents" = "Fail - Both Parents", - "other" = "Other") - ) + - ggplot2::labs( - title = "Trio Mendelian Error Distribution", - subtitle = paste0("Trios with Genotype Data Tested: ", n_total, - "\n \n", threshold_label), - x = "Mendelian Error (%)", - y = "Number of Trios", - fill = "Status" - ) + - ggplot2::theme_classic(base_size = 13) + - ggplot2::theme(legend.position = "top") - print(p) - } - } - - #### Compile and return named list #### - output_list <- base::list( - pass = final_df[status == "pass"], - fail = final_df[status == "fail"], - low_markers = final_df[status == "low_markers"], - no_genotype_data = final_df[status == "no_genotype_data"], - founders = final_df[status == "founders"], - missing_parents = final_df[status %in% c("missing_both_parents", - "missing_male_parent", - "missing_female_parent")], - full_results = final_df, - corrected_pedigree = corrected_pedigree, - plot = p - ) - return(base::invisible(output_list)) -} diff --git a/README.md b/README.md index e913a9c..a7c1735 100644 --- a/README.md +++ b/README.md @@ -21,8 +21,13 @@ BIGr is an R package developed by [Breeding Insight](https://breedinginsight.org ## Installation -To install BIGr, you'll need to have `BiocManager` installed. +The stable version of BIGr is now available on CRAN. To install from the R terminal: +```R +install.packages("BIGr") +``` + +To install the development version of BIGr, you'll need to have `BiocManager` installed. Then, install from GitHub. ```R if (!require("BiocManager", quietly = TRUE)) install.packages("BiocManager") @@ -31,9 +36,11 @@ if (!require("BiocManager", quietly = TRUE)) BiocManager::install("Breeding-Insight/BIGr", dependencies = TRUE) library(BIGr) ``` +> Note: This GitHub version of BIGr is in development. So, there could be bugs present, and the stable version of BIGr on CRAN should be viewed as more reliable. + ## Funding -BIGr development is supported by [Breeding Insight](https://breedinginsight.org/), a USDA-funded initiative based at Cornell University. +BIGr development is supported by [Breeding Insight](https://breedinginsight.org/), a USDA-funded initiative based at the University of Florida - IFAS. ## Citation diff --git a/dev/dev_history.R b/dev/dev_history.R index 34e9324..ec7f969 100644 --- a/dev/dev_history.R +++ b/dev/dev_history.R @@ -2,7 +2,7 @@ # Update dependencies in DESCRIPTION # install.packages('attachment', repos = 'https://thinkr-open.r-universe.dev') -#attachment::att_amend_desc() +attachment::att_amend_desc() # Check package coverage covr::package_coverage() diff --git a/man/allele_freq_poly.Rd b/man/allele_freq_poly.Rd deleted file mode 100644 index 6d56f79..0000000 --- a/man/allele_freq_poly.Rd +++ /dev/null @@ -1,52 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/breedtools_functions.R -\name{allele_freq_poly} -\alias{allele_freq_poly} -\title{Compute allele frequencies for populations} -\usage{ -allele_freq_poly(geno, populations, ploidy = 2) -} -\arguments{ -\item{geno}{Numeric matrix of genotypes coded as dosage of allele B -\code{{0, 1, 2, ..., ploidy}}, with individuals in rows (named) and -SNPs in columns (named).} - -\item{populations}{Named list of populations, each containing a character -vector of individual IDs belonging to that population.} - -\item{ploidy}{Integer. Ploidy level of the species. Default is \code{2}.} -} -\value{ -A data frame of allele frequencies with SNPs as rows and populations -as columns. -} -\description{ -Computes allele frequencies for specified populations from SNP array data -coded as dosage of allele B. -} -\examples{ -geno_matrix <- matrix( - c(4, 1, 4, 0, - 2, 2, 1, 3, - 0, 4, 0, 4, - 3, 3, 2, 2, - 1, 4, 2, 3), - nrow = 4, ncol = 5, byrow = FALSE, - dimnames = list(paste0("Ind", 1:4), paste0("S", 1:5)) -) -pop_list <- list( - PopA = c("Ind1", "Ind2"), - PopB = c("Ind3", "Ind4") -) -allele_freqs <- allele_freq_poly(geno = geno_matrix, populations = pop_list, ploidy = 4) -print(allele_freqs) - -} -\references{ -Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP. -Estimation of genome-wide and locus-specific breed composition in pigs. -\emph{Transl Anim Sci.} 2017;1(1):36–44. -} -\author{ -Josué Chinchilla-Vargas -} diff --git a/man/check_ped.Rd b/man/check_ped.Rd deleted file mode 100644 index 85d7ac3..0000000 --- a/man/check_ped.Rd +++ /dev/null @@ -1,62 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/check_ped.R -\name{check_ped} -\alias{check_ped} -\title{Check and Correct Common Pedigree Errors} -\usage{ -check_ped( - ped.file, - seed = NULL, - verbose = TRUE, - correct_conflicting_trios = TRUE, - correct_inconsistent_sex_roles = TRUE -) -} -\arguments{ -\item{ped.file}{Path to the pedigree text file (TSV/CSV/TXT), OR a -data.frame / data.table with columns: id, male_parent, female_parent.} - -\item{seed}{Optional integer seed for reproducibility.} - -\item{verbose}{Logical. If TRUE (default), prints the report to the console.} - -\item{correct_conflicting_trios}{Logical. If TRUE (default), sets conflicting -male_parent and female_parent to 0 and collapses to one row per ID.} - -\item{correct_inconsistent_sex_roles}{Logical. If TRUE (default), sets -male_parent and female_parent to 0 for rows involving IDs found as both, -then removes any resulting exact duplicates.} -} -\value{ -An invisible named list of data frames: -\describe{ -\item{exact_duplicates}{Exact duplicate rows found in the input.} -\item{conflicting_trios}{IDs with conflicting male_parent or female_parent assignments.} -\item{inconsistent_sex_roles}{Rows where a conflicting ID appears as male_parent or female_parent.} -\item{missing_parents}{Parent IDs absent from id, added as founders.} -\item{dependencies}{Cycles detected in the pedigree. Must be resolved manually.} -\item{corrected_pedigree}{Corrected pedigree table.} -} -} -\description{ -Reads a 3-column pedigree file (id, male_parent, female_parent) and performs -quality checks, optionally correcting detected errors. Exact duplicates and -missing parents are always corrected. Conflicting trios and inconsistent sex -roles are corrected when their respective arguments are TRUE. Cycles are -reported only and must be resolved manually. -} -\examples{ -ped_file <- system.file("check_ped_test.txt", package = "BIGr") -ped_errors <- check_ped(ped.file = ped_file, seed = 101919, verbose = FALSE) - -# Also accepts a data.table directly -library(data.table) -ped_dt <- data.table(id = c("A","B","C"), - male_parent = c("0","0","A"), - female_parent = c("0","0","B")) -ped_errors <- check_ped(ped.file = ped_dt, verbose = FALSE) - -} -\author{ -Josue Chinchilla-Vargas -} diff --git a/man/find_parentage.Rd b/man/find_parentage.Rd deleted file mode 100644 index 533572f..0000000 --- a/man/find_parentage.Rd +++ /dev/null @@ -1,76 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/find_parentage.R -\name{find_parentage} -\alias{find_parentage} -\title{Find Parentage Assignments for Progeny} -\usage{ -find_parentage( - genotypes_file, - parents_file, - progeny_file, - method = "best_pair", - min_markers = 10, - error_threshold = 5, - show_ties = TRUE, - allow_parent_selfing = FALSE, - exclude_self_match = TRUE, - verbose = TRUE, - plot_results = TRUE -) -} -\arguments{ -\item{genotypes_file}{Path to a TSV/CSV/TXT file, OR a data.frame / -data.table with an 'id' column followed by marker columns coded as 0, 1, 2.} - -\item{parents_file}{Path to a TSV/CSV/TXT file, OR a data.frame / -data.table with an 'id' column and an optional 'sex' column -('M', 'F', or 'A'). If absent, all parents are treated as ambiguous.} - -\item{progeny_file}{Path to a TSV/CSV/TXT file, OR a data.frame / -data.table with an 'id' column.} - -\item{method}{Character. One of \code{"best_male_parent"}, -\code{"best_female_parent"}, \code{"best_match"}, or -\code{"best_pair"} (default).} - -\item{min_markers}{Integer. Minimum markers required; fewer flags -\code{low_markers} (default: \code{10}).} - -\item{error_threshold}{Numeric. Maximum mismatch percentage; exceeded values -flag \code{high_error} (default: \code{5.0}). Must be between 0 and 100.} - -\item{show_ties}{Logical. If \code{TRUE}, tied best pairs are appended as -suffix columns. Default is \code{TRUE}.} - -\item{allow_parent_selfing}{Logical. If \code{FALSE}, candidate pairs with -identical male and female parent IDs are excluded. Applies only when -\code{method = "best_pair"}. Default is \code{FALSE}.} - -\item{exclude_self_match}{Logical. If \code{TRUE}, each progeny ID is -excluded from its own candidate parent set, preventing self-matches when -progeny are also present in the parents file. Default is \code{TRUE}.} - -\item{verbose}{Logical. If \code{TRUE}, prints progress and summary. -Default is \code{TRUE}.} - -\item{plot_results}{Logical. If \code{TRUE}, plots the Mendelian error -distribution. Requires \code{ggplot2}. Default is \code{TRUE}.} -} -\value{ -A named list (returned invisibly) with elements: -\describe{ -\item{pass}{Progeny with a confident parentage assignment.} -\item{high_error}{Progeny whose best assignment exceeds the error threshold.} -\item{low_markers}{Progeny with insufficient markers for a valid assignment.} -\item{full_results}{Complete data.table with all progeny and all output columns.} -\item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} -} -} -\description{ -Assigns the most likely parent(s) to each progeny from SNP genotype data -using Mendelian error rates or homozygous mismatch rates. Parents or progeny -absent from the genotype file are removed with a warning. -} -\author{ -Josue Chinchilla-Vargas -} diff --git a/man/solve_composition_poly.Rd b/man/solve_composition_poly.Rd deleted file mode 100644 index 128f7ad..0000000 --- a/man/solve_composition_poly.Rd +++ /dev/null @@ -1,81 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/breedtools_functions.R -\name{solve_composition_poly} -\alias{solve_composition_poly} -\title{Compute genome-wide breed composition} -\usage{ -solve_composition_poly( - Y, - X, - ped = NULL, - groups = NULL, - mia = FALSE, - sire = FALSE, - dam = FALSE, - ploidy = 2 -) -} -\arguments{ -\item{Y}{Numeric matrix of genotypes with individuals as rows and SNPs as -columns, coded as dosage of allele B \code{{0, 1, 2, ..., ploidy}}.} - -\item{X}{Numeric matrix of allele frequencies with SNPs as rows and -breeds/populations as columns.} - -\item{ped}{Optional data frame with pedigree information formatted with -columns \code{ID}, \code{Sire}, and \code{Dam}. When supplied, \code{QPsolve_par} is used -and only animals with genotyped parents are included.} - -\item{groups}{Optional named list of IDs grouped by breed/population. -When supplied, results are returned as a named list partitioned by group.} - -\item{mia}{Logical. Only applies when \code{ped} is supplied. If \code{TRUE}, returns -the inferred maternally inherited allele per locus per animal. Default \code{FALSE}.} - -\item{sire}{Logical. Only applies when \code{ped} is supplied. If \code{TRUE}, returns -sire genotypes per locus per animal. Default \code{FALSE}.} - -\item{dam}{Logical. Only applies when \code{ped} is supplied. If \code{TRUE}, returns -dam genotypes per locus per animal. Default \code{FALSE}.} - -\item{ploidy}{Integer. Ploidy level of the species. Default is \code{2}.} -} -\value{ -A data frame, or a named list of data frames when \code{groups} is -supplied, containing breed/ancestry composition estimates. -} -\description{ -Estimates genome-wide breed/ancestry composition for a batch of animals -using quadratic programming, with optional pedigree-assisted and -grouped-output modes. -} -\examples{ -allele_freqs_matrix <- matrix( - c(0.625, 0.500, - 0.500, 0.500, - 0.500, 0.500, - 0.750, 0.500, - 0.625, 0.625), - nrow = 5, ncol = 2, byrow = TRUE, - dimnames = list(paste0("SNP", 1:5), c("VarA", "VarB")) -) -val_geno_matrix <- matrix( - c(2, 1, 2, 3, 4, - 3, 4, 2, 3, 0), - nrow = 2, ncol = 5, byrow = TRUE, - dimnames = list(paste0("Test", 1:2), paste0("SNP", 1:5)) -) -composition <- solve_composition_poly(Y = val_geno_matrix, - X = allele_freqs_matrix, - ploidy = 4) -print(composition) - -} -\references{ -Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP. -Estimation of genome-wide and locus-specific breed composition in pigs. -\emph{Transl Anim Sci.} 2017;1(1):36–44. -} -\author{ -Josué Chinchilla-Vargas -} diff --git a/man/validate_pedigree.Rd b/man/validate_pedigree.Rd deleted file mode 100644 index 3daa9d1..0000000 --- a/man/validate_pedigree.Rd +++ /dev/null @@ -1,70 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/validate_pedigree.R -\name{validate_pedigree} -\alias{validate_pedigree} -\title{Validate Pedigree Trios Using Mendelian Error Analysis} -\usage{ -validate_pedigree( - pedigree_file, - genotypes_file, - founders_file = NULL, - trio_error_threshold = 5, - min_markers = 10, - single_parent_error_threshold = 2, - verbose = TRUE, - plot_results = TRUE -) -} -\arguments{ -\item{pedigree_file}{Path to the pedigree file (TSV/CSV/TXT), OR a -data.frame / data.table with columns: id, male_parent, female_parent.} - -\item{genotypes_file}{Path to the genotypes file (TSV/CSV/TXT), OR a -data.frame / data.table with an id column followed by marker columns -coded as 0, 1, 2.} - -\item{founders_file}{Character, optional. Path to a one-column file listing -founder IDs. Founders with both parents coded as 0 are left unchanged. -Defaults to NULL.} - -\item{trio_error_threshold}{Numeric. Maximum Mendelian error percentage to -classify a trio as pass (default: 5.0). Must be between 0 and 100.} - -\item{min_markers}{Integer. Minimum non-missing markers required to evaluate -a trio (default: 10).} - -\item{single_parent_error_threshold}{Numeric. Maximum homozygous-marker -mismatch percentage for a parent to be considered acceptable -(default: 2.0). Must be between 0 and 100.} - -\item{verbose}{Logical. If TRUE, prints progress, summary, and results to -the console (default: TRUE).} - -\item{plot_results}{Logical. If TRUE, prints a histogram of trio Mendelian -error percentages with a threshold line (default: TRUE).} -} -\value{ -An invisible named list with the following elements: -\describe{ -\item{pass}{Trios that passed the Mendelian error threshold.} -\item{fail}{Trios that failed the Mendelian error threshold.} -\item{low_markers}{Trios with insufficient markers for evaluation.} -\item{no_genotype_data}{Trios absent from the genotype file.} -\item{founders}{Trios identified as founders.} -\item{missing_parents}{Trios with one or both parents coded as 0 (non-founders).} -\item{full_results}{Complete data.table with all trios and all output columns.} -\item{corrected_pedigree}{Pedigree table after applying recommended corrections.} -\item{plot}{ggplot object if plot_results = TRUE, otherwise NULL.} -} -} -\description{ -Validates parent-offspring trios against SNP genotype data using Mendelian -error rates. Identifies incorrect parentage assignments, suggests -best-matching replacements, and outputs a corrected pedigree. Founder trios -(both parents coded as 0) are preserved unchanged if a founders file is -supplied. Trios absent from the genotype file are retained as -no_genotype_data. -} -\author{ -Josue Chinchilla-Vargas -} diff --git a/tests/testthat/test-breedtools_poly.R b/tests/testthat/test-breedtools_poly.R deleted file mode 100644 index 91ed356..0000000 --- a/tests/testthat/test-breedtools_poly.R +++ /dev/null @@ -1,31 +0,0 @@ -context("BreedTools") - - -test_that("test breedtools poly",{ - #Input variables - ref_file <- system.file("test_ref.txt", package="BIGr") - val_file <- system.file("test_test.txt", package="BIGr") - ref_ids <- system.file("ref_ids.txt", package="BIGr") - - #import files - reference = read.table(ref_file, header = T, row.names = 1, sep = "\t") - validation = read.table(val_file, header = T, row.names = 1, sep = "\t") - reference_ids = read.table(ref_ids, header = T, sep = "\t") - - #Calculations - ref_ids = lapply(as.list(reference_ids),as.character) - - freq = allele_freq_poly(reference, ref_ids, ploidy = 4) - - prediction = as.data.frame(solve_composition_poly(validation,freq, ploidy = 4)) - - #Check - freq_mean <- round(mean(as.numeric(freq)),6) - pred_mean <- round(mean(as.numeric(prediction$R2)),6) - - - expect_equal(freq_mean, 0.888889, tolerance = 0.01) - expect_equal(pred_mean, 0.841454, tolerance = 0.01) - expect_true(nrow(prediction) == 175) - -}) diff --git a/tests/testthat/test-check_ped.R b/tests/testthat/test-check_ped.R deleted file mode 100644 index ff150c3..0000000 --- a/tests/testthat/test-check_ped.R +++ /dev/null @@ -1,388 +0,0 @@ -# tests/testthat/test-check_ped.R -library(testthat) - -write_ped <- function(df) { - f <- tempfile(fileext = ".txt") - utils::write.table(df, f, sep = "\t", row.names = FALSE, quote = FALSE) - f -} - -context("check_ped – Pedigree Quality Checks") - -test_that("check_ped returns a named list of length 6", { - ped <- data.frame( - id = c("A", "B", "C"), - male_parent = c("0", "A", "A"), - female_parent = c("0", "0", "0") - ) - out <- check_ped(write_ped(ped), seed = 1, verbose = FALSE) - expect_type(out, "list") - expect_length(out, 6) - expect_named(out, c( - "exact_duplicates", - "conflicting_trios", - "inconsistent_sex_roles", - "missing_parents", - "dependencies", - "corrected_pedigree" - )) -}) - -test_that("check_ped report components are data.frames", { - ped <- data.frame( - id = c("A", "B", "C"), - male_parent = c("0", "A", "A"), - female_parent = c("0", "0", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_true(is.data.frame(out$exact_duplicates)) - expect_true(is.data.frame(out$conflicting_trios)) - expect_true(is.data.frame(out$inconsistent_sex_roles)) - expect_true(is.data.frame(out$missing_parents)) - expect_true(is.data.frame(out$dependencies)) - expect_true(is.data.frame(out$corrected_pedigree)) -}) - -test_that("corrected_pedigree has lowercase column names and no row_number", { - ped <- data.frame( - id = c("A", "B"), - male_parent = c("0", "A"), - female_parent = c("0", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_true(all(c("id", "male_parent", "female_parent") %in% - names(out$corrected_pedigree))) - expect_false("row_number" %in% names(out$corrected_pedigree)) -}) - -test_that("clean pedigree produces no issues", { - ped <- data.frame( - id = c("G1", "G2", "P1"), - male_parent = c("0", "0", "G1"), - female_parent = c("0", "0", "G2") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_equal(nrow(out$exact_duplicates), 0) - expect_equal(nrow(out$conflicting_trios), 0) - expect_equal(nrow(out$inconsistent_sex_roles), 0) - expect_equal(nrow(out$missing_parents), 0) - expect_equal(nrow(out$dependencies), 0) - expect_equal(nrow(out$corrected_pedigree), 3) -}) - -test_that("check_ped detects exact duplicates", { - ped <- data.frame( - id = c("A", "A", "B"), - male_parent = c("0", "0", "A"), - female_parent = c("0", "0", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_equal(nrow(out$exact_duplicates), 2) - expect_true(all(out$exact_duplicates$id == "A")) -}) - -test_that("exact duplicates are collapsed to one row in corrected_pedigree", { - ped <- data.frame( - id = c("A", "A", "B"), - male_parent = c("0", "0", "A"), - female_parent = c("0", "0", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_equal(sum(out$corrected_pedigree$id == "A"), 1) -}) - -test_that("check_ped detects conflicting trios", { - ped <- data.frame( - id = c("A", "A", "B"), - male_parent = c("X", "Y", "A"), - female_parent = c("M", "M", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_equal(nrow(out$conflicting_trios), 2) - expect_true(all(out$conflicting_trios$id == "A")) -}) - -test_that("correct_conflicting_trios = TRUE: conflicting field -> '0', consistent kept", { - ped <- data.frame( - id = c("A", "A", "B"), - male_parent = c("X", "Y", "A"), - female_parent = c("M", "M", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE, - correct_conflicting_trios = TRUE) - a_row <- out$corrected_pedigree[out$corrected_pedigree$id == "A", ] - expect_equal(nrow(a_row), 1) - expect_equal(a_row$male_parent, "0") - expect_equal(a_row$female_parent, "M") -}) - -test_that("correct_conflicting_trios = FALSE leaves conflicting rows as-is", { - ped <- data.frame( - id = c("A", "A", "B"), - male_parent = c("X", "Y", "A"), - female_parent = c("M", "M", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE, - correct_conflicting_trios = FALSE) - expect_equal(sum(out$corrected_pedigree$id == "A"), 2) -}) - -test_that("check_ped detects missing parents", { - ped <- data.frame( - id = c("A", "B"), - male_parent = c("0", "X"), - female_parent = c("0", "Y") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_equal(nrow(out$missing_parents), 2) - expect_true("X" %in% out$missing_parents$id) - expect_true("Y" %in% out$missing_parents$id) -}) - -test_that("missing parents are added as founder rows in corrected_pedigree", { - ped <- data.frame( - id = c("A", "B"), - male_parent = c("0", "X"), - female_parent = c("0", "Y") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_true("X" %in% out$corrected_pedigree$id) - expect_true("Y" %in% out$corrected_pedigree$id) - x_row <- out$corrected_pedigree[out$corrected_pedigree$id == "X", ] - y_row <- out$corrected_pedigree[out$corrected_pedigree$id == "Y", ] - expect_equal(x_row$male_parent, "0") - expect_equal(x_row$female_parent, "0") - expect_equal(y_row$male_parent, "0") - expect_equal(y_row$female_parent, "0") -}) - -test_that("individual that is its own parent is logged as a dependency", { - ped <- data.frame( - id = c("A", "B"), - male_parent = c("A", "0"), - female_parent = c("0", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_gt(nrow(out$dependencies), 0) -}) - -# inconsistent_sex_roles returns rows where the conflicting ID appears -# AS A PARENT — check male_parent and female_parent columns [1] -test_that("check_ped detects inconsistent sex roles", { - ped <- data.frame( - id = c("child1", "child2", "P"), - male_parent = c("P", "0", "0"), - female_parent = c("0", "P", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_true("inconsistent_sex_roles" %in% names(out)) - expect_gt(nrow(out$inconsistent_sex_roles), 0) - expect_true(any(out$inconsistent_sex_roles$male_parent == "P" | - out$inconsistent_sex_roles$female_parent == "P")) -}) - -test_that("correct_inconsistent_sex_roles = TRUE zeros out conflicting parent references", { - ped <- data.frame( - id = c("child1", "child2", "P"), - male_parent = c("P", "0", "0"), - female_parent = c("0", "P", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE, - correct_inconsistent_sex_roles = TRUE) - corr <- out$corrected_pedigree - expect_false(any(corr$male_parent == "P")) - expect_false(any(corr$female_parent == "P")) -}) - -test_that("correct_inconsistent_sex_roles = FALSE leaves conflicting references", { - ped <- data.frame( - id = c("child1", "child2", "P"), - male_parent = c("P", "0", "0"), - female_parent = c("0", "P", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE, - correct_inconsistent_sex_roles = FALSE) - corr <- out$corrected_pedigree - expect_true(any(corr$male_parent == "P" | corr$female_parent == "P")) -}) - -test_that("check_ped detects a direct two-node cycle", { - ped <- data.frame( - id = c("A", "B"), - male_parent = c("B", "A"), - female_parent = c("0", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_gt(nrow(out$dependencies), 0) -}) - -test_that("cycle-involved IDs are still present in corrected_pedigree", { - ped <- data.frame( - id = c("A", "B"), - male_parent = c("B", "A"), - female_parent = c("0", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_true("A" %in% out$corrected_pedigree$id) - expect_true("B" %in% out$corrected_pedigree$id) -}) - -test_that("check_ped errors when required columns are missing", { - bad_df <- data.frame( - animal_id = c("a", "b"), - parent1 = c("0", "a"), - parent2 = c("0", "0") - ) - expect_error( - check_ped(write_ped(bad_df), verbose = FALSE), - regexp = "missing required column" - ) -}) - -test_that("check_ped accepts mixed-case column names (ID, Male_Parent, Female_Parent)", { - ped <- data.frame( - ID = c("A", "B", "C"), - Male_Parent = c("0", "A", "A"), - Female_Parent = c("0", "0", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_length(out, 6) - expect_true(all(c("id", "male_parent", "female_parent") %in% - names(out$corrected_pedigree))) -}) - -test_that("check_ped accepts all-uppercase column names", { - ped <- data.frame( - ID = c("A", "B"), - MALE_PARENT = c("0", "A"), - FEMALE_PARENT = c("0", "0") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_length(out, 6) -}) - -test_that("check_ped accepts columns in any order", { - ped <- data.frame( - female_parent = c("0", "0"), - male_parent = c("0", "A"), - id = c("A", "B") - ) - out <- check_ped(write_ped(ped), verbose = FALSE) - expect_length(out, 6) - expect_equal(nrow(out$corrected_pedigree), 2) -}) - -test_that("seed produces reproducible results", { - ped <- data.frame( - id = c("A", "B", "C"), - male_parent = c("0", "A", "A"), - female_parent = c("0", "0", "0") - ) - f <- write_ped(ped) - out1 <- check_ped(f, seed = 42, verbose = FALSE) - out2 <- check_ped(f, seed = 42, verbose = FALSE) - expect_identical(out1$corrected_pedigree, out2$corrected_pedigree) -}) - -test_that("verbose = FALSE suppresses console output", { - ped <- data.frame( - id = c("A", "B"), - male_parent = c("0", "A"), - female_parent = c("0", "0") - ) - expect_silent(check_ped(write_ped(ped), verbose = FALSE)) -}) - -test_that("check_ped returns invisibly", { - ped <- data.frame( - id = c("A", "B"), - male_parent = c("0", "A"), - female_parent = c("0", "0") - ) - expect_invisible(check_ped(write_ped(ped), verbose = FALSE)) -}) - -test_that("no output files are written to disk", { - tmp_dir <- tempfile() - dir.create(tmp_dir) - old_wd <- getwd() - setwd(tmp_dir) - on.exit({ setwd(old_wd); unlink(tmp_dir, recursive = TRUE) }, add = TRUE) - ped <- data.frame( - id = c("A", "B"), - male_parent = c("0", "A"), - female_parent = c("0", "0") - ) - check_ped(write_ped(ped), verbose = FALSE) - expect_length(list.files(tmp_dir), 0) -}) - -# ============================================================================== -# Integration test -# Fixture has sire/dam columns renamed to male_parent/female_parent [1] -# janitor::clean_names() handles any remaining capitalization variants -# ============================================================================== -test_that("integration test with bundled fixture file", { - ped_file <- system.file("check_ped_test.txt", package = "BIGr") - skip_if(ped_file == "", "Bundled fixture file not found; skipping.") - - out <- check_ped(ped_file, seed = 101919, verbose = FALSE) - - expect_length(out, 6) - expect_gt(nrow(out$inconsistent_sex_roles), 0) - - # inconsistent_sex_roles stores rows where the conflicting ID appears - # AS A PARENT in male_parent or female_parent columns [1] - conflicting_ids <- unique(c( - out$inconsistent_sex_roles$male_parent, - out$inconsistent_sex_roles$female_parent - )) - expect_true(any(c("grandfather2", "grandfather3") %in% conflicting_ids)) - expect_equal(nrow(out$missing_parents), 8) -}) -# ============================================================================== -# In-memory input — data.frame / data.table accepted directly -# ============================================================================== - -test_that("check_ped accepts a data.frame directly", { - ped <- data.frame( - id = c("A", "B", "C"), - male_parent = c("0", "A", "A"), - female_parent = c("0", "0", "0") - ) - out <- check_ped(ped, verbose = FALSE) - expect_length(out, 6) - expect_true(all(c("id", "male_parent", "female_parent") %in% - names(out$corrected_pedigree))) -}) - -test_that("check_ped accepts a data.table directly", { - ped <- data.table::data.table( - id = c("A", "B", "C"), - male_parent = c("0", "A", "A"), - female_parent = c("0", "0", "0") - ) - out <- check_ped(ped, verbose = FALSE) - expect_length(out, 6) - expect_true(all(c("id", "male_parent", "female_parent") %in% - names(out$corrected_pedigree))) -}) - -test_that("in-memory and file-path inputs produce identical corrected_pedigree", { - ped <- data.frame( - id = c("A", "B", "C"), - male_parent = c("0", "A", "A"), - female_parent = c("0", "0", "0") - ) - out_file <- check_ped(write_ped(ped), verbose = FALSE) - out_mem <- check_ped(ped, verbose = FALSE) - expect_identical(out_file$corrected_pedigree, - out_mem$corrected_pedigree) -}) - -test_that("invalid input type raises a descriptive error for check_ped", { - expect_error( - check_ped(list(id = "A"), verbose = FALSE), - regexp = "file path" - ) -}) diff --git a/tests/testthat/test-find_parentage.R b/tests/testthat/test-find_parentage.R deleted file mode 100644 index b2c384d..0000000 --- a/tests/testthat/test-find_parentage.R +++ /dev/null @@ -1,762 +0,0 @@ -# tests/testthat/test-find_parentage.R -# Run with: testthat::test_file("tests/testthat/test-find_parentage.R") - -library(testthat) -library(data.table) -set.seed(101919) -# ============================================================================== -# Helpers -# ============================================================================== - -make_files <- function(genos, parents, progeny, dir = tempdir()) { - geno_file <- file.path(dir, paste0("genos_", sample(1e6, 1), ".txt")) - parent_file <- file.path(dir, paste0("parents_", sample(1e6, 1), ".txt")) - progeny_file <- file.path(dir, paste0("progeny_", sample(1e6, 1), ".txt")) - data.table::fwrite(genos, geno_file, sep = "\t") - data.table::fwrite(parents, parent_file, sep = "\t") - data.table::fwrite(progeny, progeny_file, sep = "\t") - list(g = geno_file, p = parent_file, pr = progeny_file) -} - -# ============================================================================== -# Base fixtures -# ============================================================================== - -# S1/D1: all 0 -> child1 (all 0) is a perfect Mendelian child of S1 x D1 -# S2/D2: all 2 -> child2 (all 2) is a perfect Mendelian child of S2 x D2 -base_genos <- data.table::data.table( - id = c("S1", "S2", "D1", "D2", "child1", "child2"), - M1 = c(0L, 2L, 0L, 2L, 0L, 2L), - M2 = c(0L, 2L, 0L, 2L, 0L, 2L), - M3 = c(0L, 2L, 0L, 2L, 0L, 2L), - M4 = c(0L, 2L, 0L, 2L, 0L, 2L), - M5 = c(0L, 2L, 0L, 2L, 0L, 2L) -) - -base_parents <- data.table::data.table(id = c("S1", "S2", "D1", "D2"), - sex = c("M", "M", "F", "F")) -child1_progeny <- data.table::data.table(id = "child1") -child2_progeny <- data.table::data.table(id = "child2") -base_progeny <- data.table::data.table(id = c("child1", "child2")) - -# Distinct marker patterns that catch row/vector recycling errors. -alignment_genos <- data.table::data.table( - id = c("M1", "F1", "M2", "F2", "C1", "C2", "C3"), - snp01 = c(0L, 0L, 2L, 2L, 0L, 2L, 1L), - snp02 = c(0L, 0L, 2L, 0L, 0L, 1L, 0L), - snp03 = c(0L, 2L, 0L, 2L, 1L, 1L, 1L), - snp04 = c(0L, 2L, 0L, 0L, 1L, 0L, 0L), - snp05 = c(2L, 0L, 1L, 1L, 1L, 1L, 1L), - snp06 = c(2L, 0L, 1L, 0L, 1L, 1L, 1L), - snp07 = c(2L, 2L, 0L, 1L, 2L, 0L, 2L), - snp08 = c(2L, 2L, 2L, 2L, 2L, 2L, 2L), - snp09 = c(0L, 1L, 0L, 2L, 1L, 1L, 1L), - snp10 = c(0L, 1L, 2L, 0L, 0L, 1L, 0L), - snp11 = c(2L, 1L, 1L, 2L, 1L, 2L, 2L), - snp12 = c(2L, 1L, 0L, 1L, 2L, 1L, 1L) -) -alignment_parents <- data.table::data.table(id = c("M1", "M2", "F1", "F2"), - sex = c("M", "M", "F", "F")) -alignment_progeny <- data.table::data.table(id = c("C1", "C2", "C3")) - -# All-zero genotypes -- every sex-compatible pair ties at 0% error. -tied_genos <- data.table::data.table( - id = c("S1", "S2", "D1", "D2", "child_tie"), - M1 = c(0L, 0L, 0L, 0L, 0L), - M2 = c(0L, 0L, 0L, 0L, 0L) -) -tied_parents <- data.table::data.table(id = c("S1", "S2", "D1", "D2"), - sex = c("M", "M", "F", "F")) -tied_progeny <- data.table::data.table(id = "child_tie") - -# Progeny are also candidate parents. -self_match_genos <- data.table::data.table( - id = c("child_self", "alt_parent"), - M1 = c(0L, 2L), - M2 = c(0L, 2L), - M3 = c(0L, 2L), - M4 = c(0L, 2L), - M5 = c(0L, 2L) -) -self_match_parents <- data.table::data.table(id = c("child_self", "alt_parent"), - sex = c("A", "A")) -self_match_progeny <- data.table::data.table(id = "child_self") - -self_pair_genos <- data.table::data.table( - id = c("child_self", "other_male", "female_zero"), - M1 = c(0L, 2L, 0L), - M2 = c(0L, 2L, 0L), - M3 = c(0L, 2L, 0L), - M4 = c(0L, 2L, 0L), - M5 = c(0L, 2L, 0L) -) -self_pair_parents <- data.table::data.table(id = c("child_self", "other_male", "female_zero"), - sex = c("M", "M", "F")) -self_pair_progeny <- data.table::data.table(id = "child_self") - -# ============================================================================== -# 1. Input validation -# ============================================================================== - -test_that("invalid method throws an error", { - f <- make_files(base_genos, base_parents, child1_progeny) - expect_error( - find_parentage(f$g, f$p, f$pr, method = "bad_method", - verbose = FALSE, plot_results = FALSE), - regexp = "Method must be one of" - ) -}) - -test_that("min_markers < 1 throws an error", { - f <- make_files(base_genos, base_parents, child1_progeny) - expect_error( - find_parentage(f$g, f$p, f$pr, min_markers = 0, - verbose = FALSE, plot_results = FALSE), - regexp = "min_markers" - ) -}) - -test_that("error_threshold out of range throws an error", { - f <- make_files(base_genos, base_parents, child1_progeny) - expect_error( - find_parentage(f$g, f$p, f$pr, error_threshold = 150, - verbose = FALSE, plot_results = FALSE), - regexp = "error_threshold" - ) - expect_error( - find_parentage(f$g, f$p, f$pr, error_threshold = -1, - verbose = FALSE, plot_results = FALSE), - regexp = "error_threshold" - ) -}) - -test_that("missing genotype file throws an error", { - f <- make_files(base_genos, base_parents, child1_progeny) - expect_error( - find_parentage("nonexistent.txt", f$p, f$pr, - verbose = FALSE, plot_results = FALSE), - regexp = "Error reading input files" - ) -}) - -test_that("parent IDs absent from genotype file raise a warning and are dropped", { - extra_parents <- rbind(base_parents, - data.table::data.table(id = "GHOST", sex = "M")) - f <- make_files(base_genos, extra_parents, child1_progeny) - expect_warning( - find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = FALSE), - regexp = "GHOST" - ) -}) - -test_that("progeny IDs absent from genotype file raise a warning and are dropped", { - extra_progeny <- rbind(child1_progeny, - data.table::data.table(id = "GHOST_KID")) - f <- make_files(base_genos, base_parents, extra_progeny) - expect_warning( - find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = FALSE), - regexp = "GHOST_KID" - ) -}) - -test_that("no valid progeny candidates after filtering stops with an error", { - ghost_progeny <- data.table::data.table(id = "NOBODY") - f <- make_files(base_genos, base_parents, ghost_progeny) - expect_warning( - expect_error( - find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = FALSE), - regexp = "No valid progeny" - ) - ) -}) - -test_that("missing sex column raises a warning and defaults to ambiguous", { - parents_no_sex <- data.table::data.table(id = c("S1", "D1")) - f <- make_files(base_genos, parents_no_sex, child1_progeny) - expect_warning( - find_parentage(f$g, f$p, f$pr, method = "best_match", - verbose = FALSE, plot_results = FALSE), - regexp = "sex" - ) -}) - -# ============================================================================== -# 2. Return structure -- named list -# ============================================================================== - -test_that("find_parentage returns an invisible named list with all required elements", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_type(out, "list") - expect_named(out, c("pass", "high_error", "low_markers", "full_results", "plot"), - ignore.order = TRUE) -}) - -test_that("full_results is a data.table", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_s3_class(out$full_results, "data.table") -}) - -test_that("full_results has one row per progeny", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(nrow(out$full_results), 1L) -}) - -test_that("named list subsets together cover all full_results rows", { - f <- make_files(base_genos, base_parents, base_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - subset_total <- nrow(out$pass) + nrow(out$high_error) + nrow(out$low_markers) - expect_equal(subset_total, nrow(out$full_results)) -}) - -test_that("plot element is NULL when plot_results = FALSE", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = FALSE) - expect_null(out$plot) -}) - -test_that("best_pair full_results has expected lowercase columns", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_true(all(c("id", "male_parent", "female_parent", - "mendelian_error_pct", "markers_tested", - "status") %in% names(out$full_results))) -}) - -test_that("best_male_parent full_results has expected lowercase columns", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_male_parent", - verbose = FALSE, plot_results = FALSE) - expect_true(all(c("id", "best_match", "mendelian_error_pct", - "markers_tested", "status") %in% names(out$full_results))) - expect_false("male_parent" %in% names(out$full_results)) -}) - -test_that("best_female_parent full_results has expected lowercase columns", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_female_parent", - verbose = FALSE, plot_results = FALSE) - expect_true(all(c("id", "best_match", "mendelian_error_pct", - "markers_tested", "status") %in% names(out$full_results))) -}) - -test_that("best_match full_results has expected lowercase columns", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_match", - verbose = FALSE, plot_results = FALSE) - expect_true(all(c("id", "best_match", "mendelian_error_pct", - "markers_tested", "status") %in% names(out$full_results))) -}) - -# ============================================================================== -# 3. Biological correctness -# ============================================================================== - -test_that("best_pair correctly identifies S1 x D1 for child1 with 0% error", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$male_parent, "S1") - expect_equal(out$full_results$female_parent, "D1") - expect_equal(as.numeric(out$full_results$mendelian_error_pct), 0) -}) - -test_that("best_pair correctly identifies S2 x D2 for child2 with 0% error", { - f <- make_files(base_genos, base_parents, child2_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$male_parent, "S2") - expect_equal(out$full_results$female_parent, "D2") - expect_equal(as.numeric(out$full_results$mendelian_error_pct), 0) -}) - -test_that("best_pair marker alignment correctly identifies all non-uniform trios", { - f <- make_files(alignment_genos, alignment_parents, alignment_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results[id == "C1"]$male_parent, "M1") - expect_equal(out$full_results[id == "C1"]$female_parent, "F1") - expect_equal(out$full_results[id == "C2"]$male_parent, "M2") - expect_equal(out$full_results[id == "C2"]$female_parent, "F2") - expect_equal(out$full_results[id == "C3"]$male_parent, "M1") - expect_equal(out$full_results[id == "C3"]$female_parent, "F2") - expect_true(all(as.numeric(out$full_results$mendelian_error_pct) == 0)) -}) - -test_that("best_male_parent identifies S1 as best male for child1", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_male_parent", - verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$best_match, "S1") -}) - -test_that("best_female_parent identifies D1 as best female for child1", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_female_parent", - verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$best_match, "D1") -}) - -test_that("single-parent marker alignment compares matching marker columns", { - f <- make_files(alignment_genos, alignment_parents, data.table::data.table(id = "C1")) - out <- find_parentage(f$g, f$p, f$pr, method = "best_match", - verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$best_match, "M1") - expect_equal(out$full_results$mendelian_error_pct, 0) -}) - -test_that("both child1 and child2 correctly assigned when run together", { - f <- make_files(base_genos, base_parents, base_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(nrow(out$full_results), 2L) - expect_equal(out$full_results[id == "child1"]$male_parent, "S1") - expect_equal(out$full_results[id == "child1"]$female_parent, "D1") - expect_equal(out$full_results[id == "child2"]$male_parent, "S2") - expect_equal(out$full_results[id == "child2"]$female_parent, "D2") -}) - -test_that("markers_tested equals number of non-NA marker columns", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$markers_tested, ncol(base_genos) - 1L) -}) - -test_that("mendelian_error_pct is between 0 and 100", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - pct <- as.numeric(out$full_results$mendelian_error_pct) - expect_true(all(pct >= 0 & pct <= 100, na.rm = TRUE)) -}) - -# ============================================================================== -# 4. status -- lowercase values -# ============================================================================== - -test_that("status = pass for perfect trio within thresholds", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - min_markers = 3, error_threshold = 5.0, - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$status, "pass") -}) - -test_that("status = low_markers when min_markers exceeds available markers", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - min_markers = 99999, error_threshold = 5.0, - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$status, "low_markers") -}) - -test_that("status = high_error when error rate exceeds threshold", { - high_error_genos <- data.table::data.table( - id = c("S1", "D1", "bad_child"), - M1 = c(0L, 0L, 2L), - M2 = c(0L, 0L, 2L), - M3 = c(0L, 0L, 2L), - M4 = c(0L, 0L, 2L), - M5 = c(0L, 0L, 2L) - ) - parents <- data.table::data.table(id = c("S1", "D1"), sex = c("M", "F")) - progeny <- data.table::data.table(id = "bad_child") - f <- make_files(high_error_genos, parents, progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - min_markers = 3, error_threshold = 5.0, - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$status, "high_error") -}) - -test_that("status column is present and lowercase in all methods", { - f <- make_files(base_genos, base_parents, child1_progeny) - for (m in c("best_pair", "best_male_parent", "best_female_parent", "best_match")) { - out <- find_parentage(f$g, f$p, f$pr, method = m, - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_true("status" %in% names(out$full_results)) - expect_true(all(out$full_results$status %in% c("pass", "high_error", "low_markers", NA))) - } -}) - -test_that("low_markers is flagged for single-parent methods too", { - f <- make_files(base_genos, base_parents, child1_progeny) - for (m in c("best_male_parent", "best_female_parent", "best_match")) { - out <- find_parentage(f$g, f$p, f$pr, method = m, - min_markers = 99999, verbose = FALSE, - plot_results = FALSE) - expect_equal(out$full_results$status, "low_markers") - } -}) - -test_that("high_error list element contains only high_error rows", { - high_error_genos <- data.table::data.table( - id = c("S1", "D1", "bad_child"), - M1 = c(0L, 0L, 2L), M2 = c(0L, 0L, 2L), M3 = c(0L, 0L, 2L), - M4 = c(0L, 0L, 2L), M5 = c(0L, 0L, 2L) - ) - parents <- data.table::data.table(id = c("S1", "D1"), sex = c("M", "F")) - progeny <- data.table::data.table(id = "bad_child") - f <- make_files(high_error_genos, parents, progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - min_markers = 3, error_threshold = 5.0, - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - if (nrow(out$high_error) > 0) - expect_true(all(out$high_error$status == "high_error")) -}) - -# ============================================================================== -# 5. Self matching controls -# ============================================================================== - -test_that("allow_parent_selfing = FALSE removes identical parent pairs from candidates", { - ambig_parents <- data.table::data.table(id = c("S1", "D1"), sex = c("A", "A")) - f <- make_files(base_genos, ambig_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - allow_parent_selfing = FALSE, show_ties = TRUE, - verbose = FALSE, plot_results = FALSE) - r <- out$full_results - if (!is.na(r$male_parent) && !is.na(r$female_parent)) - expect_false(r$male_parent == r$female_parent) -}) - -test_that("allow_parent_selfing = TRUE allows identical parent pairs", { - one_parent_genos <- data.table::data.table( - id = c("P1", "child_selfed"), - M1 = c(0L, 0L), M2 = c(0L, 0L), M3 = c(0L, 0L) - ) - one_parent <- data.table::data.table(id = "P1", sex = "A") - progeny <- data.table::data.table(id = "child_selfed") - f <- make_files(one_parent_genos, one_parent, progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - allow_parent_selfing = TRUE, - verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$male_parent, "P1") - expect_equal(out$full_results$female_parent, "P1") -}) - -test_that("default allow_parent_selfing = FALSE stops when only self-pairs are possible", { - one_parent_genos <- data.table::data.table( - id = c("P1", "child_selfed"), - M1 = c(0L, 0L), M2 = c(0L, 0L), M3 = c(0L, 0L) - ) - one_parent <- data.table::data.table(id = "P1", sex = "A") - progeny <- data.table::data.table(id = "child_selfed") - f <- make_files(one_parent_genos, one_parent, progeny) - expect_error( - find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = FALSE), - regexp = "No valid parent pairs" - ) -}) - -test_that("exclude_self_match = TRUE prevents best_match from assigning progeny to itself", { - f <- make_files(self_match_genos, self_match_parents, self_match_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_match", - exclude_self_match = TRUE, - min_markers = 3, - verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$best_match, "alt_parent") - expect_equal(out$full_results$status, "high_error") -}) - -test_that("exclude_self_match = FALSE permits best_match to assign progeny to itself", { - f <- make_files(self_match_genos, self_match_parents, self_match_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_match", - exclude_self_match = FALSE, - verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$best_match, "child_self") - expect_equal(out$full_results$mendelian_error_pct, 0) -}) - -test_that("exclude_self_match = TRUE removes parent pairs containing the progeny", { - f <- make_files(self_pair_genos, self_pair_parents, self_pair_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - exclude_self_match = TRUE, - min_markers = 3, - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_false(out$full_results$male_parent == "child_self") - expect_false(out$full_results$female_parent == "child_self") - expect_equal(out$full_results$male_parent, "other_male") - expect_equal(out$full_results$female_parent, "female_zero") - expect_equal(out$full_results$status, "high_error") -}) - -test_that("exclude_self_match = FALSE permits parent pairs containing the progeny", { - f <- make_files(self_pair_genos, self_pair_parents, self_pair_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - exclude_self_match = FALSE, - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$male_parent, "child_self") - expect_equal(out$full_results$female_parent, "female_zero") - expect_equal(as.numeric(out$full_results$mendelian_error_pct), 0) -}) - -test_that("duplicate parent rows do not create duplicate tie columns", { - duplicate_parents <- data.table::data.table(id = c("S1", "S1", "D1", "D1"), - sex = c("M", "M", "F", "F")) - f <- make_files(base_genos, duplicate_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = TRUE, verbose = FALSE, plot_results = FALSE) - expect_false(any(grepl("^male_parent_\\d", names(out$full_results)))) - expect_false(any(grepl("^female_parent_\\d", names(out$full_results)))) -}) - -# ============================================================================== -# 6. show_ties -# ============================================================================== - -test_that("show_ties = TRUE produces lowercase suffixed columns when ties exist", { - f <- make_files(tied_genos, tied_parents, tied_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = TRUE, verbose = FALSE, plot_results = FALSE) - expect_true(any(grepl("^male_parent_", names(out$full_results)))) - expect_true(any(grepl("^female_parent_", names(out$full_results)))) -}) - -test_that("show_ties = FALSE warns about ties and returns single-result columns", { - f <- make_files(tied_genos, tied_parents, tied_progeny) - expect_warning( - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE), - regexp = "tied" - ) - expect_true("male_parent" %in% names(out$full_results)) - expect_true("female_parent" %in% names(out$full_results)) - expect_false(any(grepl("^male_parent_\\d", names(out$full_results)))) - expect_false(any(grepl("^female_parent_\\d", names(out$full_results)))) -}) - -test_that("base columns are always populated even when ties exist", { - f <- make_files(tied_genos, tied_parents, tied_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = TRUE, verbose = FALSE, plot_results = FALSE) - expect_false(is.na(out$full_results$male_parent[1])) - expect_false(is.na(out$full_results$female_parent[1])) -}) - -# ============================================================================== -# 7. verbose -# ============================================================================== - -test_that("verbose = TRUE returns the result as a named list", { - f <- make_files(base_genos, base_parents, child1_progeny) - invisible(capture.output( - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = TRUE, plot_results = FALSE) - )) - expect_type(out, "list") - expect_named(out, c("pass", "high_error", "low_markers", "full_results", "plot"), - ignore.order = TRUE) -}) - -test_that("verbose = FALSE suppresses console output", { - f <- make_files(base_genos, base_parents, child1_progeny) - expect_silent( - find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = FALSE) - ) -}) - -# ============================================================================== -# 8. No write logic -- find_parentage does not write files -# ============================================================================== - -test_that("no output files are written to disk", { - tmp_dir <- tempfile() - dir.create(tmp_dir) - old_wd <- getwd() - setwd(tmp_dir) - on.exit({ setwd(old_wd); unlink(tmp_dir, recursive = TRUE) }, add = TRUE) - - f <- make_files(base_genos, base_parents, child1_progeny, dir = tmp_dir) - find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = FALSE) - - written_files <- list.files(tmp_dir, pattern = "\\.txt$|\\.jpg$|\\.csv$") - expect_equal(length(written_files), 3L) -}) - -# ============================================================================== -# 9. Sex-based candidate filtering -# ============================================================================== - -test_that("best_male_parent only assigns M or A parents", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_male_parent", - verbose = FALSE, plot_results = FALSE) - valid_male_parents <- base_parents[sex %in% c("M", "A")]$id - expect_true(out$full_results$best_match %in% valid_male_parents) -}) - -test_that("best_female_parent only assigns F or A parents", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_female_parent", - verbose = FALSE, plot_results = FALSE) - valid_female_parents <- base_parents[sex %in% c("F", "A")]$id - expect_true(out$full_results$best_match %in% valid_female_parents) -}) - -test_that("best_pair male slot contains only M or A parents", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - valid_males <- base_parents[sex %in% c("M", "A")]$id - expect_true(out$full_results$male_parent %in% valid_males) -}) - -test_that("best_pair female slot contains only F or A parents", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - valid_females <- base_parents[sex %in% c("F", "A")]$id - expect_true(out$full_results$female_parent %in% valid_females) -}) - -# ============================================================================== -# 10. Edge cases -# ============================================================================== - -test_that("single progeny individual is handled correctly", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_equal(nrow(out$full_results), 1L) -}) - -test_that("all-NA marker column does not cause an error", { - na_genos <- data.table::copy(base_genos) - na_genos[, M1 := NA_integer_] - f <- make_files(na_genos, base_parents, child1_progeny) - expect_no_error( - find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = FALSE) - ) -}) - -test_that("id column in full_results contains the correct progeny IDs", { - f <- make_files(base_genos, base_parents, child1_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_setequal(out$full_results$id, child1_progeny$id) -}) - -test_that("multiple progeny are all represented in full_results", { - f <- make_files(base_genos, base_parents, base_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - expect_setequal(out$full_results$id, base_progeny$id) -}) - -test_that("single parent pair does not cause dimension errors", { - single_pair_parents <- data.table::data.table(id = c("S1", "D1"), - sex = c("M", "F")) - f <- make_files(base_genos, single_pair_parents, child1_progeny) - expect_no_error( - find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, plot_results = FALSE) - ) -}) - -test_that("one row returned per progeny for single-parent methods", { - f <- make_files(base_genos, base_parents, child1_progeny) - for (m in c("best_male_parent", "best_female_parent", "best_match")) { - out <- find_parentage(f$g, f$p, f$pr, method = m, - verbose = FALSE, plot_results = FALSE) - expect_equal(nrow(out$full_results), 1L) - } -}) - -test_that("all candidate parents excluded by exclude_self_match yields low_markers", { - only_self_parents <- data.table::data.table(id = "child_self", sex = "A") - f <- make_files(self_match_genos, only_self_parents, self_match_progeny) - out <- find_parentage(f$g, f$p, f$pr, method = "best_match", - exclude_self_match = TRUE, - verbose = FALSE, plot_results = FALSE) - expect_true(is.na(out$full_results$best_match)) - expect_equal(out$full_results$markers_tested, 0L) - expect_equal(out$full_results$status, "low_markers") -}) - -# ============================================================================== -# 11. plot element -# ============================================================================== - -test_that("plot element is a ggplot when plot_results = TRUE", { - skip_if_not_installed("ggplot2") - f <- make_files(base_genos, base_parents, child1_progeny) - out <- suppressWarnings( - find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = TRUE) - ) - expect_s3_class(out$plot, "ggplot") -}) - -# ============================================================================== -# 12. Return value is invisible -# ============================================================================== - -test_that("find_parentage returns invisibly", { - f <- make_files(base_genos, base_parents, child1_progeny) - expect_invisible( - find_parentage(f$g, f$p, f$pr, method = "best_pair", - verbose = FALSE, plot_results = FALSE) - ) -}) - -# ============================================================================== -# 13. In-memory input -- data.frame / data.table accepted directly -# ============================================================================== - -test_that("find_parentage accepts data.table inputs directly", { - out <- find_parentage(base_genos, base_parents, child1_progeny, - method = "best_pair", show_ties = FALSE, - verbose = FALSE, plot_results = FALSE) - expect_s3_class(out$full_results, "data.table") - expect_equal(nrow(out$full_results), 1L) -}) - -test_that("find_parentage accepts data.frame inputs directly", { - out <- find_parentage(as.data.frame(base_genos), - as.data.frame(base_parents), - as.data.frame(child1_progeny), - method = "best_pair", show_ties = FALSE, - verbose = FALSE, plot_results = FALSE) - expect_s3_class(out$full_results, "data.table") - expect_equal(nrow(out$full_results), 1L) -}) - -test_that("in-memory and file-path inputs produce identical results for find_parentage", { - f <- make_files(base_genos, base_parents, child1_progeny) - out_file <- find_parentage(f$g, f$p, f$pr, method = "best_pair", - show_ties = FALSE, verbose = FALSE, - plot_results = FALSE) - out_mem <- find_parentage(base_genos, base_parents, child1_progeny, - method = "best_pair", show_ties = FALSE, - verbose = FALSE, plot_results = FALSE) - expect_equal(out_file$full_results$male_parent, - out_mem$full_results$male_parent) - expect_equal(out_file$full_results$mendelian_error_pct, - out_mem$full_results$mendelian_error_pct) - expect_equal(out_file$full_results$status, - out_mem$full_results$status) -}) - -test_that("invalid input type raises a descriptive error for find_parentage", { - expect_error( - find_parentage(list(id = "S1"), base_parents, child1_progeny, - verbose = FALSE, plot_results = FALSE), - regexp = "Error reading input files" - ) -}) diff --git a/tests/testthat/test-validate_pedigree.R b/tests/testthat/test-validate_pedigree.R deleted file mode 100644 index 7930ae2..0000000 --- a/tests/testthat/test-validate_pedigree.R +++ /dev/null @@ -1,749 +0,0 @@ -# tests/testthat/test-validate_pedigree.R -# Run with: testthat::test_file("tests/testthat/test-validate_pedigree.R") - -library(testthat) -library(data.table) - -# ============================================================================== -# Helpers -# ============================================================================== - -make_genos <- function() { - n_markers <- 20 - marker_names <- paste0("M", seq_len(n_markers)) - dt <- data.table( - id = c("IND_A", "IND_B", "IND_C", "IND_D", "IND_E"), - rbind( - rep(0L, n_markers), # IND_A — all ref homozygous - rep(2L, n_markers), # IND_B — all alt homozygous - rep(1L, n_markers), # IND_C — all het: valid child of IND_A x IND_B - rep(0L, n_markers), # IND_D — impossible child of IND_B x IND_A - rep(0L, n_markers) # IND_E — all ref - ) - ) - setnames(dt, c("id", marker_names)) - dt -} - -make_pedigree <- function() { - # IND_C: perfect Mendelian child of IND_A x IND_B -> pass - # IND_D: declared parents swapped -> fail - data.table( - id = c("IND_C", "IND_D"), - male_parent = c("IND_A", "IND_B"), - female_parent = c("IND_B", "IND_A") - ) -} - -write_temp_files <- function(genos = make_genos(), ped = make_pedigree()) { - ped_file <- tempfile(fileext = ".txt") - genos_file <- tempfile(fileext = ".txt") - fwrite(ped, ped_file, sep = "\t") - fwrite(genos, genos_file, sep = "\t") - list(ped = ped_file, genos = genos_file) -} - -# ============================================================================== -# 1. Input validation -# ============================================================================== - -test_that("trio_error_threshold out of range raises an error", { - f <- write_temp_files() - expect_error( - validate_pedigree(f$ped, f$genos, trio_error_threshold = 150, - verbose = FALSE, plot_results = FALSE), - regexp = "trio_error_threshold" - ) - expect_error( - validate_pedigree(f$ped, f$genos, trio_error_threshold = -1, - verbose = FALSE, plot_results = FALSE), - regexp = "trio_error_threshold" - ) -}) - -test_that("single_parent_error_threshold out of range raises an error", { - f <- write_temp_files() - expect_error( - validate_pedigree(f$ped, f$genos, single_parent_error_threshold = 101, - verbose = FALSE, plot_results = FALSE), - regexp = "single_parent_error_threshold" - ) - expect_error( - validate_pedigree(f$ped, f$genos, single_parent_error_threshold = -5, - verbose = FALSE, plot_results = FALSE), - regexp = "single_parent_error_threshold" - ) -}) - -test_that("boundary values 0 and 100 are accepted for trio_error_threshold", { - f <- write_temp_files() - expect_no_error( - validate_pedigree(f$ped, f$genos, trio_error_threshold = 0, - verbose = FALSE, plot_results = FALSE) - ) - expect_no_error( - validate_pedigree(f$ped, f$genos, trio_error_threshold = 100, - verbose = FALSE, plot_results = FALSE) - ) -}) - -test_that("nonexistent pedigree file throws 'Error reading input files'", { - f <- write_temp_files() - expect_error( - validate_pedigree("nonexistent.txt", f$genos, - verbose = FALSE, plot_results = FALSE), - regexp = "Error reading input files" - ) -}) - -test_that("nonexistent genotypes file throws 'Error reading input files'", { - f <- write_temp_files() - expect_error( - validate_pedigree(f$ped, "nonexistent.txt", - verbose = FALSE, plot_results = FALSE), - regexp = "Error reading input files" - ) -}) - -test_that("missing required pedigree column raises an error", { - bad_ped <- data.table(id = "IND_C", parent1 = "IND_A", female_parent = "IND_B") - f <- write_temp_files(ped = bad_ped) - expect_error( - validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE), - regexp = "missing required columns" - ) -}) - -test_that("missing id column in genotypes raises an error", { - bad_genos <- copy(make_genos()) - setnames(bad_genos, "id", "SampleID") - f <- write_temp_files(genos = bad_genos) - expect_error( - validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE), - regexp = "id" - ) -}) - -test_that("all trios with no genotype data stops with an error", { - ped <- data.table(id = "GHOST", male_parent = "IND_A", female_parent = "IND_B") - f <- write_temp_files(ped = ped) - expect_error( - validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE), - regexp = "No valid trios remain" - ) -}) - -test_that("unreadable founders file raises an error", { - f <- write_temp_files() - expect_error( - validate_pedigree(f$ped, f$genos, founders_file = "nonexistent_founders.txt", - verbose = FALSE, plot_results = FALSE) - ) -}) - -# ============================================================================== -# 2. Return structure — named list -# ============================================================================== - -test_that("returns an invisible named list with all required elements", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_type(out, "list") - expect_named(out, c("pass", "fail", "low_markers", "no_genotype_data", - "founders", "missing_parents", "full_results", - "corrected_pedigree", "plot"), - ignore.order = TRUE) -}) - -test_that("validate_pedigree returns invisibly", { - f <- write_temp_files() - expect_invisible( - validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - ) -}) - -test_that("full_results is a data.table", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_s3_class(out$full_results, "data.table") -}) - -test_that("full_results has one row per pedigree entry", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_equal(nrow(out$full_results), 2L) -}) - -test_that("full_results has all expected lowercase columns", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expected_cols <- c( - "id", "orig_male_parent", "orig_female_parent", - "trio_mendelian_error_pct", "trio_markers_tested", "status", - "recommended_correction", - "male_parent_hom_error_pct", "female_parent_hom_error_pct", - "best_male_candidate", "best_male_candidate_error_pct", - "best_female_candidate", "best_female_candidate_error_pct" - ) - expect_true(all(expected_cols %in% names(out$full_results))) -}) - -test_that("corrected_pedigree is a data.table with lowercase columns", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_s3_class(out$corrected_pedigree, "data.table") - expect_true(all(c("id", "male_parent", "female_parent") %in% - names(out$corrected_pedigree))) -}) - -test_that("corrected_pedigree has same number of rows as original pedigree", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_equal(nrow(out$corrected_pedigree), nrow(make_pedigree())) -}) - -test_that("plot element is NULL when plot_results = FALSE", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_null(out$plot) -}) - -test_that("plot element is a ggplot when plot_results = TRUE", { - skip_if_not_installed("ggplot2") - f <- write_temp_files() - out <- suppressWarnings( - validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = TRUE) - ) - expect_s3_class(out$plot, "ggplot") -}) - -test_that("named list subsets sum correctly to full_results row count", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - subset_total <- nrow(out$pass) + nrow(out$fail) + nrow(out$low_markers) + - nrow(out$no_genotype_data) + nrow(out$founders) + nrow(out$missing_parents) - expect_equal(subset_total, nrow(out$full_results)) -}) - -# ============================================================================== -# 3. pass / fail / low_markers statuses -# ============================================================================== - -test_that("pass trio is correctly identified with 0% error", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_C"] - expect_equal(nrow(r), 1L) - expect_equal(r$status, "pass") - expect_equal(r$trio_mendelian_error_pct, 0) - expect_equal(r$recommended_correction, "none") -}) - -test_that("fail trio is correctly identified with error above threshold", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_D"] - expect_equal(nrow(r), 1L) - expect_equal(r$status, "fail") - expect_gt(r$trio_mendelian_error_pct, 5.0) -}) - -test_that("fail trio has a non-NA recommended_correction", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_D"] - expect_false(is.na(r$recommended_correction)) - expect_false(r$recommended_correction == "none") -}) - -test_that("fail trio with one acceptable parent gets remove_* correction", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_D"] - expect_true(r$recommended_correction %in% - c("remove_male_parent", "remove_female_parent", "remove_both", - "keep_both")) -}) - -test_that("trio_mendelian_error_pct is 0 for a perfect Mendelian trio", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results[id == "IND_C"]$trio_mendelian_error_pct, 0) -}) - -test_that("trio_mendelian_error_pct is between 0 and 100 for all trios", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - pct <- out$full_results$trio_mendelian_error_pct - expect_true(all(pct >= 0 & pct <= 100, na.rm = TRUE)) -}) - -test_that("trio_markers_tested equals number of markers for complete data", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results[id == "IND_C"]$trio_markers_tested, 20L) -}) - -test_that("low_markers status assigned when markers_tested < min_markers", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, - plot_results = FALSE, min_markers = 25L) - expect_true(all(out$full_results$status == "low_markers")) - expect_true(all(grepl("^low_markers_", out$full_results$recommended_correction))) -}) - -test_that("NA markers reduce trio_markers_tested and do not cause errors", { - genos <- make_genos() - genos[id == "IND_C", M1 := NA_integer_] - genos[id == "IND_C", M2 := NA_integer_] - f <- write_temp_files(genos = genos) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results[id == "IND_C"]$trio_markers_tested, 18L) - expect_equal(out$full_results[id == "IND_C"]$status, "pass") -}) - -test_that("pass list element contains only pass rows", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - if (nrow(out$pass) > 0) - expect_true(all(out$pass$status == "pass")) -}) - -test_that("fail list element contains only fail rows", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - if (nrow(out$fail) > 0) - expect_true(all(out$fail$status == "fail")) -}) - -test_that("low_markers list element contains only low_markers rows", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, - plot_results = FALSE, min_markers = 25L) - if (nrow(out$low_markers) > 0) - expect_true(all(out$low_markers$status == "low_markers")) -}) - -test_that("raising trio_error_threshold turns fail rows into pass rows", { - f <- write_temp_files() - strict <- validate_pedigree(f$ped, f$genos, trio_error_threshold = 5.0, - verbose = FALSE, plot_results = FALSE) - lenient <- validate_pedigree(f$ped, f$genos, trio_error_threshold = 100.0, - verbose = FALSE, plot_results = FALSE) - expect_gte(nrow(lenient$pass), nrow(strict$pass)) -}) - -# ============================================================================== -# 4. missing parent statuses -# ============================================================================== - -test_that("missing_male_parent status and recommendation are correct", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_E", male_parent = "0", - female_parent = "IND_B")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_E"] - expect_equal(r$status, "missing_male_parent") - expect_equal(r$recommended_correction, "none") - expect_false(is.na(r$best_male_candidate)) - expect_true(is.na(r$best_female_candidate)) -}) - -test_that("missing_female_parent status and recommendation are correct", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_E", male_parent = "IND_A", - female_parent = "0")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_E"] - expect_equal(r$status, "missing_female_parent") - expect_equal(r$recommended_correction, "none") - expect_true(is.na(r$best_male_candidate)) - expect_false(is.na(r$best_female_candidate)) -}) - -test_that("missing_both_parents status and recommendations are correct", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_E", male_parent = "0", - female_parent = "0")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_E"] - expect_equal(r$status, "missing_both_parents") - expect_equal(r$recommended_correction, "none") - expect_false(is.na(r$best_male_candidate)) - expect_false(is.na(r$best_female_candidate)) -}) - -test_that("best_male_candidate for missing_male_parent excludes the known female parent", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_E", male_parent = "0", - female_parent = "IND_B")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_E"] - expect_false(r$best_male_candidate == "IND_B") -}) - -test_that("best_female_candidate for missing_female_parent excludes the known male parent", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_E", male_parent = "IND_A", - female_parent = "0")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_E"] - expect_false(r$best_female_candidate == "IND_A") -}) - -test_that("missing_parents list element contains only missing_* rows", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_E", male_parent = "0", - female_parent = "0")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - if (nrow(out$missing_parents) > 0) - expect_true(all(grepl("^missing_", out$missing_parents$status))) -}) - -# ============================================================================== -# 5. founders status -# ============================================================================== - -test_that("founders status is assigned when ID is in founders list with 0 0 parents", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_A", male_parent = "0", - female_parent = "0")) - founders_file <- tempfile(fileext = ".txt") - fwrite(data.table(id = "IND_A"), founders_file, - sep = "\t", quote = FALSE, col.names = FALSE) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, founders_file = founders_file, - verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_A"] - expect_equal(r$status, "founders") - expect_equal(r$recommended_correction, "none") - expect_true(is.na(r$best_male_candidate)) - expect_true(is.na(r$best_female_candidate)) -}) - -test_that("founders list element contains only founders rows", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_A", male_parent = "0", - female_parent = "0")) - founders_file <- tempfile(fileext = ".txt") - fwrite(data.table(id = "IND_A"), founders_file, - sep = "\t", quote = FALSE, col.names = FALSE) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, founders_file = founders_file, - verbose = FALSE, plot_results = FALSE) - if (nrow(out$founders) > 0) - expect_true(all(out$founders$status == "founders")) -}) - -test_that("non-founder rows still evaluated normally when founders file is supplied", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_A", male_parent = "0", - female_parent = "0")) - founders_file <- tempfile(fileext = ".txt") - fwrite(data.table(id = "IND_A"), founders_file, - sep = "\t", quote = FALSE, col.names = FALSE) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, founders_file = founders_file, - verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results[id == "IND_C"]$status, "pass") - expect_equal(out$full_results[id == "IND_D"]$status, "fail") -}) - -test_that("0 0 parents NOT in founders list get missing_both_parents", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_E", male_parent = "0", - female_parent = "0")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results[id == "IND_E"]$status, "missing_both_parents") -}) - -test_that("founder row does not appear in pass, fail, or missing_parents", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_A", male_parent = "0", - female_parent = "0")) - founders_file <- tempfile(fileext = ".txt") - fwrite(data.table(id = "IND_A"), founders_file, - sep = "\t", quote = FALSE, col.names = FALSE) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, founders_file = founders_file, - verbose = FALSE, plot_results = FALSE) - expect_false("IND_A" %in% out$pass$id) - expect_false("IND_A" %in% out$fail$id) - expect_false("IND_A" %in% out$missing_parents$id) -}) - -# ============================================================================== -# 6. no_genotype_data status -# ============================================================================== - -test_that("no_genotype_data is flagged for progeny absent from genotype file", { - ped <- rbind(make_pedigree(), - data.table(id = "GHOST", male_parent = "IND_A", - female_parent = "IND_B")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "GHOST"] - expect_equal(nrow(r), 1L) - expect_equal(r$status, "no_genotype_data") - expect_equal(r$recommended_correction, "none") -}) - -test_that("no_genotype_data rows have NA/0 for all analysis columns", { - ped <- rbind(make_pedigree(), - data.table(id = "GHOST", male_parent = "IND_A", - female_parent = "IND_B")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "GHOST"] - expect_true(is.na(r$trio_mendelian_error_pct)) - expect_equal(r$trio_markers_tested, 0L) - expect_true(is.na(r$best_male_candidate)) - expect_true(is.na(r$best_female_candidate)) -}) - -test_that("no_genotype_data flagged when a declared parent is absent from genotype file", { - ped <- rbind(make_pedigree(), - data.table(id = "IND_C_GHOST", male_parent = "GHOST_DAD", - female_parent = "IND_B")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results[id == "IND_C_GHOST"]$status, "no_genotype_data") -}) - -test_that("no_genotype_data list element contains only no_genotype_data rows", { - ped <- rbind(make_pedigree(), - data.table(id = "GHOST", male_parent = "IND_A", - female_parent = "IND_B")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - if (nrow(out$no_genotype_data) > 0) - expect_true(all(out$no_genotype_data$status == "no_genotype_data")) -}) - -test_that("valid trios still evaluated correctly when ghost rows are present", { - ped <- rbind(make_pedigree(), - data.table(id = "GHOST", male_parent = "IND_A", - female_parent = "IND_B")) - f <- write_temp_files(ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results[id == "IND_C"]$status, "pass") - expect_equal(out$full_results[id == "IND_D"]$status, "fail") -}) - -# ============================================================================== -# 7. corrected_pedigree contents -# ============================================================================== - -test_that("corrected_pedigree: pass parents are unchanged", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - corr <- out$corrected_pedigree - expect_equal(as.character(corr[id == "IND_C"]$male_parent), "IND_A") - expect_equal(as.character(corr[id == "IND_C"]$female_parent), "IND_B") -}) - -test_that("corrected_pedigree: removed parent set to 0 for remove_male_parent", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - r <- out$full_results[id == "IND_D"] - corr <- out$corrected_pedigree - if (r$recommended_correction == "remove_male_parent") { - expect_equal(corr[id == "IND_D"]$male_parent, "0") - } -}) - -test_that("corrected_pedigree: removed parent set to 0 for remove_female_parent", { - # construct a trio where IND_A (all 0) is correct male and female is wrong - genos <- make_genos() - ped <- data.table(id = "IND_E", male_parent = "IND_A", - female_parent = "IND_B") - # IND_E is all ref (0); IND_A is all ref (0); IND_B is all alt (2) - # IND_E as child of IND_A x IND_B is impossible → remove_female_parent - f <- write_temp_files(genos = genos, ped = ped) - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - corr <- out$corrected_pedigree - r <- out$full_results[id == "IND_E"] - if (r$recommended_correction == "remove_female_parent") { - expect_equal(corr[id == "IND_E"]$female_parent, "0") - expect_false(corr[id == "IND_E"]$male_parent == "0") - } -}) - -test_that("corrected_pedigree preserves id column values", { - f <- write_temp_files() - out <- validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_setequal(out$corrected_pedigree$id, make_pedigree()$id) -}) - -# ============================================================================== -# 8. No write logic — function does not write files -# ============================================================================== - -test_that("no output files are written to disk", { - f <- write_temp_files() - tmp_dir <- tempfile() - dir.create(tmp_dir) - old_wd <- getwd() - setwd(tmp_dir) - on.exit({ setwd(old_wd); unlink(tmp_dir, recursive = TRUE) }, add = TRUE) - validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - written_files <- list.files(tmp_dir) - expect_length(written_files, 0) -}) - -# ============================================================================== -# 9. verbose -# ============================================================================== - -test_that("verbose = FALSE suppresses console output", { - f <- write_temp_files() - expect_silent( - validate_pedigree(f$ped, f$genos, verbose = FALSE, plot_results = FALSE) - ) -}) - -test_that("verbose = TRUE returns valid named list without error", { - f <- write_temp_files() - invisible(capture.output( - out <- validate_pedigree(f$ped, f$genos, verbose = TRUE, plot_results = FALSE) - )) - expect_type(out, "list") - expect_named(out, c("pass", "fail", "low_markers", "no_genotype_data", - "founders", "missing_parents", "full_results", - "corrected_pedigree", "plot"), - ignore.order = TRUE) -}) - -# ============================================================================== -# 10. Mendelian error correctness -# ============================================================================== - -test_that("0x0 parents produce 0% error for dosage-0 child", { - genos <- data.table( - id = c("S", "D", "C"), - M1 = c(0L, 0L, 0L), M2 = c(0L, 0L, 0L), M3 = c(0L, 0L, 0L), - M4 = c(0L, 0L, 0L), M5 = c(0L, 0L, 0L) - ) - ped <- data.table(id = "C", male_parent = "S", female_parent = "D") - f <- write_temp_files(genos = genos, ped = ped) - out <- validate_pedigree(f$ped, f$genos, trio_error_threshold = 5.0, - min_markers = 1, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$trio_mendelian_error_pct, 0) -}) - -test_that("2x2 parents produce 0% error for dosage-2 child", { - genos <- data.table( - id = c("S", "D", "C"), - M1 = c(2L, 2L, 2L), M2 = c(2L, 2L, 2L), M3 = c(2L, 2L, 2L), - M4 = c(2L, 2L, 2L), M5 = c(2L, 2L, 2L) - ) - ped <- data.table(id = "C", male_parent = "S", female_parent = "D") - f <- write_temp_files(genos = genos, ped = ped) - out <- validate_pedigree(f$ped, f$genos, trio_error_threshold = 5.0, - min_markers = 1, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$trio_mendelian_error_pct, 0) -}) - -test_that("0x0 parents produce 100% error for dosage-2 child", { - genos <- data.table( - id = c("S", "D", "C"), - M1 = c(0L, 0L, 2L), M2 = c(0L, 0L, 2L), M3 = c(0L, 0L, 2L), - M4 = c(0L, 0L, 2L), M5 = c(0L, 0L, 2L) - ) - ped <- data.table(id = "C", male_parent = "S", female_parent = "D") - f <- write_temp_files(genos = genos, ped = ped) - out <- validate_pedigree(f$ped, f$genos, trio_error_threshold = 5.0, - min_markers = 1, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$trio_mendelian_error_pct, 100) -}) - -test_that("0x2 parents produce 0% error for dosage-1 child", { - genos <- data.table( - id = c("S", "D", "C"), - M1 = c(0L, 2L, 1L), M2 = c(0L, 2L, 1L), M3 = c(0L, 2L, 1L), - M4 = c(0L, 2L, 1L), M5 = c(0L, 2L, 1L) - ) - ped <- data.table(id = "C", male_parent = "S", female_parent = "D") - f <- write_temp_files(genos = genos, ped = ped) - out <- validate_pedigree(f$ped, f$genos, trio_error_threshold = 5.0, - min_markers = 1, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$trio_mendelian_error_pct, 0) -}) - -test_that("0x2 parents produce 100% error for dosage-0 child", { - genos <- data.table( - id = c("S", "D", "C"), - M1 = c(0L, 2L, 0L), M2 = c(0L, 2L, 0L), M3 = c(0L, 2L, 0L), - M4 = c(0L, 2L, 0L), M5 = c(0L, 2L, 0L) - ) - ped <- data.table(id = "C", male_parent = "S", female_parent = "D") - f <- write_temp_files(genos = genos, ped = ped) - out <- validate_pedigree(f$ped, f$genos, trio_error_threshold = 5.0, - min_markers = 1, verbose = FALSE, plot_results = FALSE) - expect_equal(out$full_results$trio_mendelian_error_pct, 100) -}) -# ============================================================================== -# 11. In-memory input — data.frame / data.table accepted directly -# ============================================================================== - -test_that("validate_pedigree accepts a data.table pedigree directly", { - f <- write_temp_files() - ped <- make_pedigree() - out <- validate_pedigree(ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_s3_class(out$full_results, "data.table") - expect_equal(nrow(out$full_results), 2L) -}) - -test_that("validate_pedigree accepts a data.table genotypes object directly", { - f <- write_temp_files() - genos <- make_genos() - out <- validate_pedigree(f$ped, genos, verbose = FALSE, plot_results = FALSE) - expect_s3_class(out$full_results, "data.table") - expect_equal(nrow(out$full_results), 2L) -}) - -test_that("validate_pedigree accepts both inputs as data.tables directly", { - ped <- make_pedigree() - genos <- make_genos() - out <- validate_pedigree(ped, genos, verbose = FALSE, plot_results = FALSE) - expect_s3_class(out$full_results, "data.table") - expect_equal(nrow(out$full_results), 2L) -}) - -test_that("validate_pedigree accepts a data.frame pedigree directly", { - f <- write_temp_files() - ped <- as.data.frame(make_pedigree()) - out <- validate_pedigree(ped, f$genos, verbose = FALSE, plot_results = FALSE) - expect_s3_class(out$full_results, "data.table") - expect_equal(nrow(out$full_results), 2L) -}) - -test_that("in-memory and file-path inputs produce identical results for validate_pedigree", { - f <- write_temp_files() - ped <- make_pedigree() - genos <- make_genos() - out_file <- validate_pedigree(f$ped, f$genos, - verbose = FALSE, plot_results = FALSE) - out_memory <- validate_pedigree(ped, genos, - verbose = FALSE, plot_results = FALSE) - expect_equal(out_file$full_results$status, - out_memory$full_results$status) - expect_equal(out_file$full_results$trio_mendelian_error_pct, - out_memory$full_results$trio_mendelian_error_pct) -}) - -test_that("invalid input type raises an error for validate_pedigree", { - f <- write_temp_files() - expect_error( - validate_pedigree(list(id = "IND_C"), f$genos, - verbose = FALSE, plot_results = FALSE), - regexp = "Error reading input files" - ) -})