diff --git a/docs/BIGpopA_user_guide.html b/docs/BIGpopA_user_guide.html new file mode 100644 index 0000000..d45ce52 --- /dev/null +++ b/docs/BIGpopA_user_guide.html @@ -0,0 +1,4503 @@ + + + + + + + + + + + + + + + +BIGpopA Tutorial + + + + + + + + + + + + + + + + + +
+ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
+ + + +
+
+
+
+
+ +
+ + + + + + + +
+

1 Introduction

+

BIGpopA is an R package for pedigree quality control, parentage +assignment and breed/line composition from SNP genotypes. It works for +diploid and polyploid species and is the analytical engine behind the Familia Shiny +app. This guide runs the functions in the order a breeding program +usually needs them:

+
    +
  1. Convert genotypes to the dosage table used by every function. -> +vcf_to_dosage(), ped_to_dosage()
    2. +
  2. Is the pedigree record clean? -> check_ped()
    3. +
  3. Do the recorded parents agree with the genotypes? -> +validate_pedigree()
    4. +
  4. Who are the parents of an individual with unknown or wrong parents? +-> find_parentage()
    5. +
  5. What is the line / breed composition of each sample? -> +allele_freq_poly(), solve_composition_poly(), +assign_composition_line()
    6. +
+ +++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
FunctionUse it toMain output
vcf_to_dosage()Convert a VCF to allele-B dosagesDosage table or matrix
ped_to_dosage()Convert a PLINK .ped (diploid) to allele-B dosagesDosage table or matrix
check_ped()Fix structural pedigree errors before any genotype workCorrected pedigree, issue tables
validate_pedigree()Test recorded trios against genotypesStatus per trio, corrected pedigree, marker summary
find_parentage()Assign or fill in parents from a candidate poolBest parent(s), error %, ties
allele_freq_poly()Allele frequencies of reference linesFrequency table (markers x lines)
solve_composition_poly()Composition of samples from reference frequenciesProportions and R2 per sample
assign_composition_line()Percentages and assigned line per sampleComposition table
+

Remarks on the data (violations give misleading results):

+
    +
  • Genotypes are allele-B dosages (0 to ploidy), not +allele pairs. Convert VCF or PLINK files with the two converters +below.
    • +
  • Ploidy is correct. The Mendelian rules depend on +the gametes each parent can produce (ploidy / 2 alleles). A wrong ploidy +produces many false failures, so set it to the ploidy of your dosages. +Even ploidies use a polysomic gamete-range test (assumes autopolyploid +inheritance, so it is conservative for allopolyploids). Odd ploidies +(e.g. triploid) fall back to an opposite-homozygote test on homozygous +markers only, which has less power.
    • +
  • Offspring come from normal sexual reproduction. +Selfed, inbred, doubled-haploid or clonal material can break the +assumptions of the trio test (see the Tips section).
    • +
  • Marker names and individual ids match exactly +between files (case and spaces matter).
    • +
+

BIGpopA does not call dosages; use your pipeline first (we recommend +BIGapp).

+
+
+

2 Install

+

BIGpopA requires R >= 4.4.0. Install it from CRAN by copying and +pasting this code into your R console:

+
install.packages("BIGpopA")
+library(BIGpopA)
+
+
+

3 Input data

+

Every function accepts either a file path or an in-memory object +(data.frame / data.table). founders_file accepts a file +path, a character vector of ids, or a data.frame whose first column +holds the ids. This guide builds a small toy dataset that exercises +every result category.

+ +++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
ObjectFormatUsed by
Pedigreeid, male_parent, +female_parent. Unknown parent = 0.check_ped(), validate_pedigree()
Genotypesid + one column per marker (allele-B dosage, +NA = missing)validate_pedigree(), find_parentage()
Parentsid and sex (M, +F, A = either). Without a sex +column all candidates are treated as A.find_parentage()
Progenyid column. For fill_pedigree: the 3-column +pedigree (unknown parent = 0, NA or +empty).find_parentage()
FoundersFile path, character vector or data.frame (first column)validate_pedigree(), find_parentage()
+
library(BIGpopA)
+
+# SNP genotypes coded 0/1/2 (NA = missing); C_low is mostly missing
+genotypes = data.frame(
+  id    = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low"),
+  snp01 = c(  0,   0,   2,   2,   0,   2,   1,     0,     NA),
+  snp02 = c(  0,   0,   2,   0,   0,   1,   0,     0,     NA),
+  snp03 = c(  0,   2,   0,   2,   1,   1,   1,     1,     NA),
+  snp04 = c(  0,   2,   0,   0,   1,   0,   0,     1,     NA),
+  snp05 = c(  2,   0,   1,   1,   1,   1,   1,     1,     NA),
+  snp06 = c(  2,   0,   1,   0,   1,   1,   1,     1,     NA),
+  snp07 = c(  2,   2,   0,   1,   2,   0,   2,     2,     NA),
+  snp08 = c(  2,   2,   2,   2,   2,   2,   2,     2,     NA),
+  snp09 = c(  0,   1,   0,   2,   1,   1,   1,     1,      0),
+  snp10 = c(  0,   1,   2,   0,   0,   1,   0,     0,      0),
+  snp11 = c(  2,   1,   1,   2,   1,   2,   2,     1,      2),
+  snp12 = c(  2,   1,   0,   1,   2,   1,   1,     2,      2),
+  stringsAsFactors = FALSE
+)
+
+# pedigree with founders coded 0, a mis-assigned trio (C_bad),
+# and a duplicated row (C_missing)
+pedigree = data.frame(
+  id            = c("M1","F1","M2","F2","C1","C2","C3","C_bad","C_low","C_missing","C_missing"),
+  male_parent   = c("0", "0", "0", "0", "M1","M2","M1","M2",   "M1",   "M1",       "M1"),
+  female_parent = c("0", "0", "0", "0", "F1","F2","F2","F1",   "F1",   "F1",       "F1"),
+  stringsAsFactors = FALSE
+)
+
+# candidate parents (with sex) and progeny to assign
+parents = data.frame(
+  id  = c("M1","M2","F1","F2"),
+  sex = c("M", "M", "F", "F"),
+  stringsAsFactors = FALSE
+)
+
+progeny = data.frame(id = c("C1","C2","C3","C_bad","C_low","C_missing"),
+                     stringsAsFactors = FALSE)
+
+pedigree
+
+ +
+

The pedigree has four founders (M1, F1, +M2, F2, both parents 0) and +several crosses, including a duplicated row (C_missing), a +mis-assigned trio (C_bad) and a sample with mostly missing +genotypes (C_low).

+
+

3.1 Converting genotype +files

+

Use these converters when your genotypes come as a VCF or a PLINK +file. Both return an id column followed by one column per +marker, ready for validate_pedigree() and +find_parentage().

+
# VCF (.vcf or .vcf.gz); ploidy is inferred from the GT field when NULL
+geno_vcf = vcf_to_dosage("my_markers.vcf.gz", ploidy = 4)
+
+# PLINK .ped (diploid only); a .map file next to it is used for marker names
+geno_ped = ped_to_dosage("my_markers.ped")
+
+# format = "matrix" returns individuals x markers, the input of
+# allele_freq_poly() and solve_composition_poly()
+geno_mat = vcf_to_dosage("my_markers.vcf.gz", ploidy = 4, format = "matrix")
+

Useful vcf_to_dosage() arguments:

+ ++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
ArgumentEffect
ploidyChecked against the ploidy observed in the file; a warning is issued +if they disagree
dosage_fieldRead a FORMAT field with dosages (e.g. "DS", +"UD") instead of counting alleles in GT; +ploidy must be given
ref_dosageCount reference instead of alternate alleles (default +FALSE)
biallelic_onlyDrop markers with more than one ALT allele (default +TRUE)
min_marker_call_rate, +min_sample_call_rateDrop markers / individuals with a lower call rate (default 0 = no +filtering)
+
+
+
+

4 Pedigree cleaning

+

check_ped() works on the columns id, +male_parent and female_parent. It applies the +checks in this order:

+ +++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
IssueMeaningAction
Exact duplicatesThe identical row appears more than onceAlways removed
Conflicting triosSame id, different parents in different rowsDisagreeing parents set to 0 and rows collapsed, or +flagged (correct_conflicting_trios = FALSE)
Missing parentsA parent is listed but has no row of its ownAlways added as a row with unknown parents (0)
Inconsistent sex rolesSame individual used as sire and as damSet to 0 as a parent everywhere, or flagged +(correct_inconsistent_sex_roles = FALSE)
Cycles / dependenciesAn individual is its own ancestor, or is recorded as its own +parentFlagged only; must be fixed manually, since there is no safe +automatic choice
+

Unknown parents must be coded 0.

+
clean_ped_results = check_ped(pedigree, verbose = FALSE)
+
+# the corrected, analysis-ready pedigree
+clean_ped = clean_ped_results$corrected_pedigree
+clean_ped
+
+ +
+

The result also holds one table per issue, for example +clean_ped_results$exact_duplicates (here the repeated +C_missing row) and +clean_ped_results$missing_parents. An empty table means the +issue was not found. Set verbose = TRUE to print a full +report to the console.

+
+
+

5 Pedigree +validation

+

validate_pedigree() compares every trio (offspring, male +parent, female parent) with the genotypes and computes a Mendelian error +rate. A marker is testable when the offspring and both parents +have a call.

+
ped_validate_results = validate_pedigree(clean_ped, genotypes,
+                                         verbose = FALSE, plot_results = TRUE)
+

+
ped_validate_report = ped_validate_results$full_results
+ped_validate_report
+
+ +
+
+

5.1 Parameters

+ +++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
ArgumentDefaultEffect
ploidy2Ploidy of the dosages
trio_error_threshold5Maximum % of testable markers with a Mendelian error for a trio to +pass
min_markers10Minimum testable markers needed to judge a trio; fewer gives +low_markers
single_parent_error_threshold2Maximum homozygous mismatch % for a single parent to pass the +single-parent check
founders_fileNULLFile path, vector of ids or data.frame; ids listed here with both +parents 0 are reported as founders
marker_trio_tableFALSEAlso return the marker x trio table (see below)
plot_resultsTRUEDraw a histogram of trio error percentages
+
+
+

5.2 Status of each +trio

+ ++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
StatusMeaning
passMendelian error at or below trio_error_threshold
failError above the threshold
low_markersFewer than min_markers testable markers
no_genotype_dataAt least one of the three individuals is not in the genotype +file
foundersBoth parents 0 and the id is in +founders_file
missing_both_parents, missing_male_parent, +missing_female_parentOne or both parents are 0
no_dataNo testable marker at all
+

In the toy data C1, C2 and C3 +pass, C_bad fails, C_low is +low_markers and C_missing has +no_genotype_data.

+
+ +
+

5.4 Marker summary

+

marker_summary lists every marker with +trios_tested, trios_mismatch and +mismatch_pct, sorted worst first. It answers a different +question than the trio table: is the problem the individual or the +marker? A few markers with a mismatch % far above the rest +(paralogous loci, poor clustering, null alleles) inflate the error of +many trios. Remove them from the genotype file and rerun.

+
head(ped_validate_results$marker_summary)
+
+ +
+
+
+

5.5 Marker x trio +table

+

Set marker_trio_table = TRUE to see, for each marker and +each trio, where an error comes from. Columns are named +id_male_female.

+
mt_results = validate_pedigree(clean_ped, genotypes, verbose = FALSE,
+                               plot_results = FALSE, marker_trio_table = TRUE)
+mt_results$marker_trio_table
+
+ +
+ ++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
LabelMeaning
matchThe trio is consistent at this marker
mismatch_male, mismatch_femaleOnly that parent is inconsistent with the offspring
mismatch_bothEach parent alone is inconsistent with the offspring
mismatch_combinationEach parent alone is compatible, but the pair cannot produce the +offspring
missing_male, missing_female, +missing_progenyThat individual has no call; every missing individual is named (for +example missing_progeny_male_female)
untestableNo call is missing, but the marker cannot be tested (odd ploidy, +offspring not homozygous)
+
+
+
+

6 Missing data

+

Missing marker calls (NA) are handled marker by marker +and never count as errors:

+
    +
  • validate_pedigree(). A marker is testable +for a trio only if the offspring and both parents have a call (odd +ploidy: the offspring is homozygous and at least one parent has a call +and is homozygous). The error % is mismatches / testable markers. Trios +with fewer testable markers than min_markers are labeled +low_markers; trios with none are no_data (only +in full_results).
    • +
  • Single-parent check. Uses only markers where both +parent and offspring are homozygous. If there are none, the parent +counts as failing.
    • +
  • Marker summary. mismatch_pct is +NA for markers not tested in any trio.
    • +
  • find_parentage(). Same testable-marker logic per +candidate; a candidate with no testable markers gets +markers_tested = 0 and low_markers.
    • +
  • allele_freq_poly(). A marker missing in every +reference individual of a line has no frequency; +solve_composition_poly() stops with an error listing these +markers. Remove them and rerun.
    • +
+
+
+

7 Parentage +assignment

+

When parentage is unknown, find_parentage() searches a +pool of candidate parents for the best assignment to each progeny.

+
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
+
+ +
+
+

7.1 Parameters

+ +++++ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +
ArgumentDefaultEffect
method"best_pair"The search (see below)
ploidy2Ploidy of the dosages
error_threshold5Error % above which the best candidate gets +high_error
min_markers10If the best candidate has fewer testable markers, the status is +low_markers
top_n1Report the best top_n candidates per progeny
show_tiesTRUEReport tied_candidates
allow_parent_selfingFALSEAllow the same individual as male and female parent
exclude_self_matchTRUEDo not let a progeny match itself as parent
founders_fileNULLFile path, vector of ids or data.frame; founders are skipped by +fill_pedigree
+ ++++ + + + + + + + + + + + + + + + + + + + + + + + + +
MethodSearches
best_pairMale x female pair with the lowest Mendelian error (needs both sexes +in the parents table)
best_male_parent, best_female_parentThe single candidate of that sex with the lowest homozygous mismatch +(compares only markers homozygous in both candidate and offspring); use +when the other parent is known or unavailable
best_matchSame homozygous mismatch test over candidates of either sex
fill_pedigreeMissing parents of a 3-column pedigree (see below)
+

The best candidate is always reported. The status is +pass (error at or below error_threshold), +high_error (above it, so treat the assignment with +caution), low_markers, or +known_parent_no_genotype. In the toy data +C_bad, which failed validation, is recovered to its true +pair (M1/F1) at 0% error, while +C_low remains low_markers.

+
+
+

7.2 Top N and ties

+

tied_candidates is the number of candidates (pairs for +best_pair) that share the lowest error of a progeny; values +above 1 mean the best assignment is ambiguous. This is common with few +markers or with clones. With top_n > 1, rank 1 stays in +the base columns and ranks 2 to top_n are added as suffix +columns (male_parent_2, female_parent_2, +mendelian_error_pct_2).

+
find_parentage(genotypes, parents, progeny, method = "best_pair", top_n = 3,
+               verbose = FALSE, plot_results = FALSE)$full_results
+
+ +
+
+
+

7.3 Filling in missing +parents of a pedigree

+

With method = "fill_pedigree", progeny is a +pedigree with id, male_parent and +female_parent (for example +validate_pedigree()$corrected_pedigree). Parents coded +0, NA or empty are unknown. A known male +parent searches only females, a known female parent only males, and rows +with both unknown get the full pair search; rows with both parents known +are skipped. The parents table (with sex M, F +or A) is still required. The output adds +search_mode and known_parent_error_pct.

+
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
+
+ +
+
+
+
+

8 Breed and line +composition

+

BIGpopA estimates genome-wide breed or line composition for diploid +and polyploid populations with a breedTools-style quadratic programming +approach (Funkhouser et al. 2017). The workflow has three steps: compute +reference-line allele frequencies, estimate the composition of each +sample, and assign each sample to a line. This example mirrors an +autotetraploid analysis (ploidy = 4) that assigns samples +to two cultivar reference panels, Jewel and Draper.

+
# 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
+
##         Jewel    Draper
+## S1 0.08333333 0.9166667
+## S2 0.08333333 0.9166667
+## S3 0.16666667 0.8333333
+## S4 0.16666667 0.8333333
+## S5 0.08333333 0.9166667
+## S6 0.16666667 0.8333333
+

The references should be representative, unadmixed individuals; the +more distinct the references, the better the estimate. Marker names must +match between the reference and the samples.

+
# samples to assign: samples x SNPs (same SNPs as the reference)
+validation = data.frame(
+  S1 = c(0, 4, 2),
+  S2 = c(1, 3, 2),
+  S3 = c(0, 4, 2),
+  S4 = c(1, 3, 2),
+  S5 = c(0, 4, 2),
+  S6 = c(1, 3, 2),
+  row.names = c("Sample1","Sample2","Sample3")
+)
+
+# composition (one column per panel plus an R2 model-fit column)
+prediction = as.data.frame(solve_composition_poly(validation, freq, ploidy = 4))
+prediction
+
+ +
+

The model returns an R2 value, but it is expected to be low because a +discrete allele dosage is regressed on a continuous allele +frequency.

+

assign_composition_line() turns the result into the +table used by Familia: the estimated percentage of each panel and a +Predicted line. A sample is assigned to its +highest-proportion line only if that proportion is at least +assign_threshold (default 50); otherwise it is labelled +Undetermined. With assign_threshold = 0 every +sample is assigned.

+
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).

+
+
+

9 Complete workflow

+
# genotypes
+genotypes = vcf_to_dosage("my_markers.vcf.gz", ploidy = 4)
+
+# pedigree quality control
+clean_ped = check_ped(pedigree)$corrected_pedigree
+validated = validate_pedigree(clean_ped, genotypes, ploidy = 4)
+assigned  = find_parentage(genotypes, parents, progeny, method = "fill_pedigree",
+                           ploidy = 4)
+
+# breed/line composition
+freq       = allele_freq_poly(reference, ref_ids, ploidy = 4)
+prediction = solve_composition_poly(validation, freq, ploidy = 4)
+composition_table = assign_composition_line(prediction, assign_threshold = 50)
+
+
+

10 Tips

+
    +
  • Typical workflow: check_ped() -> +validate_pedigree() -> find_parentage() +(fill_pedigree) on the corrected pedigree.
    • +
  • Make sure you are using the right ploidy. A wrong ploidy gives many +false failures.
    • +
  • Selfed or inbred offspring can fail the trio test; check the +single-parent errors and the unresolved flag.
    • +
  • Ties are common with few markers or with clones; always check +tied_candidates.
    • +
  • Check marker_summary before blaming individuals; remove +bad markers and rerun.
    • +
+
+
+

11 How to cite

+
+

11.1 BIGpopA

+

Chinchilla-Vargas J, Sandercock AM (2026). BIGpopA: Pedigree +Validation Genetic Composition of Diploids & Polyploids. R +package version 2.2.0. CRAN: https://CRAN.R-project.org/package=BIGpopA

+

Run citation("BIGpopA") in R for the BibTeX entry.

+
+
+

11.2 Methods

+
    +
  • Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP (2017). +Estimation of genome-wide and locus-specific breed composition in pigs. +Translational Animal Science 1(1):36-44. https://doi.org/10.2527/tas2016.0003
    • +
+
+
+
+

12 Session +information

+
sessionInfo()
+
## R version 4.5.2 (2025-10-31)
+## Platform: aarch64-apple-darwin20
+## Running under: macOS 27.0.1
+## 
+## Matrix products: default
+## BLAS:   /System/Library/Frameworks/Accelerate.framework/Versions/A/Frameworks/vecLib.framework/Versions/A/libBLAS.dylib 
+## LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1
+## 
+## locale:
+## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
+## 
+## time zone: America/Indiana/Indianapolis
+## tzcode source: internal
+## 
+## attached base packages:
+## [1] stats     graphics  grDevices utils     datasets  methods   base     
+## 
+## other attached packages:
+## [1] BIGpopA_2.2.0
+## 
+## loaded via a namespace (and not attached):
+##  [1] gtable_0.3.6        jsonlite_2.0.0      dplyr_1.2.1        
+##  [4] compiler_4.5.2      tidyselect_1.2.1    stringr_1.6.0      
+##  [7] snakecase_0.11.1    jquerylib_0.1.4     scales_1.4.0       
+## [10] yaml_2.3.12         fastmap_1.2.0       ggplot2_4.0.3      
+## [13] R6_2.6.1            labeling_0.4.3      generics_0.1.4     
+## [16] knitr_1.52          tibble_3.3.1        janitor_2.2.1      
+## [19] lubridate_1.9.5     bslib_0.12.0        pillar_1.11.1      
+## [22] RColorBrewer_1.1-3  rlang_1.3.0         stringi_1.8.9      
+## [25] cachem_1.1.0        xfun_0.61           quadprog_1.5-8     
+## [28] sass_0.4.10         S7_0.2.2            otel_0.2.0         
+## [31] timechange_0.4.0    cli_3.6.6           withr_3.0.3        
+## [34] magrittr_2.0.5      digest_0.6.39       grid_4.5.2         
+## [37] rstudioapi_0.19.0   lifecycle_1.0.5     vctrs_0.7.3        
+## [40] evaluate_1.0.5      glue_1.8.1          data.table_1.18.6.1
+## [43] farver_2.1.2        rmarkdown_2.32      tools_4.5.2        
+## [46] pkgconfig_2.0.3     htmltools_0.5.9
+
+ + + +
+
+ +
+ + + + + + + + + + + + + + + + +