--- title: "BIGpopA Tutorial" subtitle: "Pedigree cleaning, validation, and parentage assignment from SNP genotype data" author: "Breeding Insight - Josue Chinchilla-Vargas" date: "`r Sys.Date()`" output: rmarkdown::html_vignette: toc: true toc_depth: 3 vignette: > %\VignetteIndexEntry{BIGpopA Tutorial} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, echo=FALSE} knitr::opts_chunk$set( echo = TRUE, message = FALSE, warning = FALSE, fig.align = "center", fig.width = 7, fig.height = 4.5 ) ``` # Introduction This tutorial demonstrates a typical pedigree workflow using the BIGpopA R package, the analytical engine behind [Familia](https://github.com/Breeding-Insight/Familia), Breeding Insight's Shiny application for pedigree validation and ancestry assessment. We run three functions in sequence: * check_ped() — detect and correct structural pedigree errors (duplicates, conflicting trios, inconsistent sex roles, missing parents). * validate_pedigree() — test parent–offspring trios for Mendelian errors against SNP genotypes, flag incorrect assignments and summarize mismatches per marker. * find_parentage() — assign the most likely parent(s) to progeny from a pool of candidate parents, or fill in the missing parents of a pedigree. It also estimates genome-wide breed/line composition for diploid and polyploid populations with allele_freq_poly() and solve_composition_poly(), shown at the end. The methods are species-agnostic and apply to a wide range of plant and animal breeding programs. ## Install * It requires R version superior to 4.4.0 Install BIGpopA from CRAN (if not already installed): ```{r, eval=FALSE} install.packages("BIGpopA") ``` # Input data BIGpopA uses a small toy dataset that exercises every result category. Each function accepts a file path or a data.frame / data.table; to keep this vignette self-contained we build the objects inline: a pedigree, a SNP genotype matrix coded 0/1/2, and candidate parent / progeny lists. ```{r} library(BIGpopA) # SNP genotypes coded 0/1/2 (NA = missing); C_low is mostly missing genotypes = data.frame( id = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low"), snp01 = c( 0, 0, 2, 2, 0, 2, 1, 0, NA), snp02 = c( 0, 0, 2, 0, 0, 1, 0, 0, NA), snp03 = c( 0, 2, 0, 2, 1, 1, 1, 1, NA), snp04 = c( 0, 2, 0, 0, 1, 0, 0, 1, NA), snp05 = c( 2, 0, 1, 1, 1, 1, 1, 1, NA), snp06 = c( 2, 0, 1, 0, 1, 1, 1, 1, NA), snp07 = c( 2, 2, 0, 1, 2, 0, 2, 2, NA), snp08 = c( 2, 2, 2, 2, 2, 2, 2, 2, NA), snp09 = c( 0, 1, 0, 2, 1, 1, 1, 1, 0), snp10 = c( 0, 1, 2, 0, 0, 1, 0, 0, 0), snp11 = c( 2, 1, 1, 2, 1, 2, 2, 1, 2), snp12 = c( 2, 1, 0, 1, 2, 1, 1, 2, 2), stringsAsFactors = FALSE ) # pedigree with founders coded 0, a mis-assigned trio (C_bad), # and a duplicated row (C_missing) pedigree = data.frame( id = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low","C_missing","C_missing"), male_parent = c("0", "0", "0", "0", "M1","M2","M1","M2", "M1", "M1", "M1"), female_parent = c("0", "0", "0", "0", "F1","F2","F2","F1", "F1", "F1", "F1"), stringsAsFactors = FALSE ) # candidate parents (with sex) and progeny to assign parents = data.frame( id = c("M1","M2","F1","F2"), sex = c("M", "M", "F", "F"), stringsAsFactors = FALSE ) progeny = data.frame( id = c("C1","C2","C3","C_bad","C_low","C_missing"), stringsAsFactors = FALSE ) pedigree genotypes ``` The pedigree has four founders (M1, F1, M2, F2, both parents coded 0) and several crosses, including a deliberately duplicated row (C_missing), a misassigned trio (C_bad), and a sample with mostly missing genotypes (C_low). # Pedigree cleaning check_ped() operates on the pedigree columns id, male_parent, female_parent. Exact duplicates and missing parents are always corrected; conflicting trios and inconsistent sex roles are corrected by default. Cycles are reported only and must be resolved manually. Set verbose = TRUE to print a full report to the console. ```{r} #check_ped clean_ped_results = check_ped(pedigree, verbose = FALSE) # the corrected, analysis-ready pedigree clean_ped = clean_ped_results$corrected_pedigree clean_ped ``` The returned list also exposes per-issue tables for review, e.g. clean_ped_results$exact_duplicates (here, the repeated C_missing row) and clean_ped_results$missing_parents. # Pedigree validation validate_pedigree() compares each trio against the genotype matrix and computes a Mendelian error rate. Trios are classified as pass, fail, low_markers, no_genotype_data, founders, or missing_parents. To fill in missing or removed parents, run find_parentage() with method = "fill_pedigree". With plot_results = TRUE the function draws a histogram of trio error percentages. ```{r} #validate_ped ped_validate_results = validate_pedigree(clean_ped, genotypes, verbose = FALSE, plot_results = TRUE) ped_validate_report = ped_validate_results$full_results ped_validate_report ``` Reading the status column: C1, C2, C3 pass; C_bad fails (its recorded male parent shows a high homozygous-mismatch rate, so the recommended correction is to remove it); C_low is flagged low_markers because too few markers were genotyped; and C_missing has no_genotype_data. Key tuning arguments are trio_error_threshold (default 5.0), min_markers (default 10), and single_parent_error_threshold (default 2.0). The result also contains marker_summary, a per-marker table with the number of trios tested, the number with a Mendelian mismatch, and the mismatch percentage, sorted worst first. It covers every trio with both parents known and genotyped, and helps find markers (for example paralogous or miscalled SNPs) that drive many failures. ```{r} #marker_summary head(ped_validate_results$marker_summary) ``` # Parentage assignment When parentage is unknown, find_parentage() searches a pool of candidate parents for the best assignment to each progeny. The method argument selects the search: best_pair (default), best_male_parent, best_female_parent, or best_match. ```{r} #find_parentage find_parentage_results = find_parentage(genotypes, parents, progeny, method = "best_pair", verbose = FALSE, plot_results = TRUE) find_parentage_report = find_parentage_results$full_results find_parentage_report ``` Each progeny is assigned the pair minimizing Mendelian error and labeled pass, high_error, or low_markers. Notice that C_bad, which failed validation against its recorded parents, is here correctly recovered to its true pair (M1/F1) at 0% error, while C_low remains low_markers. Useful arguments include error_threshold, min_markers, show_ties, allow_parent_selfing, and exclude_self_match. The tied_candidates column gives the number of candidates (pairs for best_pair) that share each progeny's lowest error; values above 1 mean the best assignment is ambiguous. Set top_n to report the best top_n candidates per progeny; rank 1 stays in the base columns and ranks 2 to top_n are added as suffix columns (for example male_parent_2, female_parent_2, mendelian_error_pct_2). ```{r} #top_n find_parentage(genotypes, parents, progeny, method = "best_pair", top_n = 3, verbose = FALSE, plot_results = FALSE)$full_results ``` ## Filling in missing parents of a pedigree With method = "fill_pedigree", progeny is a pedigree with id, male_parent and female_parent columns (for example validate_pedigree()$corrected_pedigree). Parents coded 0, NA or empty are unknown. A known male parent searches only females, a known female parent only males, and rows with both unknown get the full pair search; rows with both parents known are skipped. The parents file (with sex M, F or A) is still required. The founders_file argument skips founders. The output adds search_mode and known_parent_error_pct. ```{r} #fill_pedigree ped_to_fill = data.frame(id = c("C1", "C2"), male_parent = c("M1", "0"), female_parent = c("0", "0")) find_parentage(genotypes, parents, ped_to_fill, method = "fill_pedigree", verbose = FALSE, plot_results = FALSE)$full_results ``` # Breed and line composition Beyond pedigree QC, BIGpopA estimates genome-wide breed or line composition for diploid and polyploid populations using a breedTools-style quadratic programming approach (Funkhouser et al. 2017). This example mirrors a blueberry (autotetraploid, ploidy = 4) analysis that assigns validation samples to two cultivar reference panels, Jewel and Draper. The workflow has two steps: compute reference-panel allele frequencies with allele_freq_poly(), then estimate each validation sample's composition with solve_composition_poly(). ```{r} library(dplyr) # reference genotypes: individuals x SNPs, tetraploid dosage of allele B (0..4) reference = data.frame( S1 = c(0, 1, 0, 4, 3, 4), S2 = c(1, 0, 0, 3, 4, 4), S3 = c(0, 1, 1, 4, 3, 3), S4 = c(1, 0, 1, 3, 4, 3), S5 = c(0, 1, 0, 4, 3, 4), S6 = c(1, 0, 1, 3, 4, 3), row.names = c("Jewel1","Jewel2","Jewel3","Draper1","Draper2","Draper3") ) # reference panel membership (which individuals define each cultivar) ref_ids = list( Jewel = c("Jewel1","Jewel2","Jewel3"), Draper = c("Draper1","Draper2","Draper3") ) # allele frequencies of the reference panels (SNPs x panels) freq = allele_freq_poly(reference, ref_ids, ploidy = 4) freq ``` ```{r} # validation samples to assign: samples x SNPs (same SNPs as the reference) validation = data.frame( S1 = c(0, 4, 2), S2 = c(1, 3, 2), S3 = c(0, 4, 2), S4 = c(1, 3, 2), S5 = c(0, 4, 2), S6 = c(1, 3, 2), row.names = c("Sample1","Sample2","Sample3") ) # breed/line composition (panel columns plus an R2 model-fit column) prediction = as.data.frame(solve_composition_poly(validation, freq, ploidy = 4)) prediction ``` It is important to note that while the model produces an R2 value, it is expected to be low since an allele dosage (discrete number) is being regressed on an allele frequency (continuous value). `assign_composition_line()` turns the result into the table used by Familia: the estimated percentage of each panel column and a `Predicted line`. A sample is assigned to its highest-proportion line only if that proportion is at least `assign_threshold` (default 50); otherwise it is labelled `Undetermined`. Use `assign_threshold = 0` to always assign the highest line. ```{r} pred_results = assign_composition_line(solve_composition_poly(validation, freq, ploidy = 4)) pred_results ``` Sample1 is assigned Jewel and Sample2 Draper (each essentially 100%), while the admixed Sample3 splits roughly evenly between them. For diploid data set `ploidy = 2`. Save the table with `write.table(pred_results, "composition.txt", sep = "\t", row.names = FALSE)`. # Complete workflow ```{r, eval=FALSE} # pedigree quality control clean_ped = check_ped(pedigree)$corrected_pedigree validated = validate_pedigree(clean_ped, genotypes)$full_results assigned = find_parentage(genotypes, parents, progeny, method = "best_pair")$full_results # breed/line composition freq = allele_freq_poly(reference, ref_ids, ploidy = 4) prediction = solve_composition_poly(validation, freq, ploidy = 4) composition_table = assign_composition_line(prediction, assign_threshold = 50) ``` # Session information ```{r} sessionInfo() ``` # How to cite If you use BIGpopA in research, run citation("BIGpopA") for the current reference. ```{r} citation("BIGpopA") ``` # Acknowledgments Developed as part of the [Breeding Insight](https://breedinginsight.org/) initiative to provide open-source, data-driven tools for modern breeding programs.