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. Read genotypes as allele-B dosages; VCF and PLINK files are converted automatically. -> vcf_to_dosage(), ped_to_dosage()
  2. Is the pedigree record clean? -> check_ped()
  3. Do the recorded parents agree with the genotypes? -> validate_pedigree()
  4. Who are the parents of an individual with unknown or wrong parents? -> find_parentage()
  5. What is the line / breed composition of each sample? -> allele_freq_poly(), solve_composition_poly(), assign_composition_line()
Function Use it to Main output
vcf_to_dosage() Convert a VCF to allele-B dosages Dosage table or matrix
ped_to_dosage() Convert a PLINK .ped (diploid) to allele-B dosages Dosage table or matrix
check_ped() Fix structural pedigree errors before any genotype work Corrected pedigree, issue tables
validate_pedigree() Test recorded trios against genotypes Status per trio, corrected pedigree, marker summary
find_parentage() Assign or fill in parents from a candidate pool Best parent(s), error %, ties
allele_freq_poly() Allele frequencies of reference lines Frequency table (markers x lines)
solve_composition_poly() Composition of samples from reference frequencies Proportions and R2 per sample
assign_composition_line() Percentages and assigned line per sample Composition table

Remarks on the data (violations give misleading results):

  • Genotypes are allele-B dosages (0 to the ploidy number), not allele pairs: the number of copies of the counted allele, so 0, 1 or 2 in a diploid and 0 to 4 in a tetraploid. VCF and PLINK files are converted automatically; vcf_to_dosage() and ped_to_dosage() show the converted table (see Converting genotype files).
  • 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 Tips).
  • Individual ids and marker names match exactly between files (case and spaces matter). Column names are more flexible; see Input data.

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

BIGpopA works with up to five inputs. Each can be a file path or an in-memory object (data.frame / data.table); genotypes can also be a VCF or a PLINK .ped file. This guide builds a small toy dataset that covers the main result categories.

Input What it holds Format Used by
Pedigree The recorded parents of each individual, one row per individual Columns id, male_parent and female_parent. These names are required but matched ignoring case, spaces and dots (ID or Male parent also work). Unknown parent = 0. check_ped(), validate_pedigree()
Genotypes The marker calls of every genotyped individual, parents and offspring in the same table id + one column per marker with the allele-B dosage (0 to ploidy), missing call = NA; or a VCF / PLINK .ped file validate_pedigree(), find_parentage()
Parents The candidates find_parentage() searches id and an optional sex (M, F, or A = either). Without sex all candidates are A. find_parentage()
Progeny The individuals whose parents you want to find id. For method = "fill_pedigree", a pedigree instead (unknown parent = 0, NA or empty). find_parentage()
Founders (optional) The base-generation individuals, whose parents are unknown by design A vector of ids, a one-column file or a data.frame (ids in the first column) validate_pedigree(), find_parentage(method = "fill_pedigree")

The same genotype table is used by validate_pedigree() and find_parentage(). Parents and founders do different jobs. parents is the pool of candidates searched by find_parentage() and can include any genotyped individual. founders only marks the individuals expected to have unknown parents: validate_pedigree() reports them as founders instead of missing_both_parents, and fill_pedigree does not search parents for them. In the toy data the four founders are also the four candidate parents.

library(BIGpopA)

# Genotypes: allele-B dosage 0/1/2 (diploid), NA = missing call.
# Parents and offspring are in the same table; 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: one row per individual, unknown parents coded 0. It has four
# founders (M1, F1, M2, F2), a mis-assigned trio (C_bad), a sample with mostly
# missing genotypes (C_low) 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
)

# Founders (optional): individuals whose parents are unknown by design
founders = c("M1", "F1", "M2", "F2")

# For find_parentage(): 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

C_bad is recorded as M2 x F1 but is really a child of M1 x F1, and C_missing is not genotyped. The next sections show how each function handles these cases.

3.1 Converting genotype files

VCF and PLINK files do not need to be converted first: validate_pedigree(), find_parentage(), allele_freq_poly() and solve_composition_poly() accept the file path and convert it with their ploidy. vcf_to_dosage() (VCF, any ploidy) and ped_to_dosage() (PLINK .ped, diploid) run the same conversion when you want to see or reuse the table. Below are three markers of the toy data for M1, F1 and C1 in each format.

example.vcf (tab-separated; the GT calls 0/0, 0/1 and 1/1 are 0, 1 and 2 copies of the ALT allele)

##fileformat=VCFv4.3
##FORMAT=<ID=GT,Number=1,Type=String,Description="Genotype">
#CHROM  POS ID  REF ALT QUAL    FILTER  INFO    FORMAT  M1  F1  C1
chr1    100 snp03   A   G   .   .   .   GT  0/0 1/1 0/1
chr1    200 snp05   C   T   .   .   .   GT  1/1 0/0 0/1
chr1    300 snp09   G   A   .   .   .   GT  0/0 0/1 0/1

example.ped (six pedigree columns, of which only the individual id in column 2 is used, then two alleles per marker)

FAM M1 0 0 1 -9 A A T T A A
FAM F1 0 0 2 -9 G G C C A G
FAM C1 M1 F1 0 -9 A G C T A G

example.map (one row per marker: chromosome, marker name, genetic position, physical position)

1 snp03 0 100
1 snp05 0 200
1 snp09 0 300

Converting either file gives the same table:

geno_vcf = vcf_to_dosage("example.vcf", verbose = FALSE)   # ploidy read from the GT calls
geno_vcf
geno_ped = ped_to_dosage("example.ped", verbose = FALSE)   # example.map in the same folder gives the marker names
geno_ped

Both have the layout of genotypes and, here, the same values:

genotypes[genotypes$id %in% c("M1", "F1", "C1"), c("id", "snp03", "snp05", "snp09")]

In a .ped file the counted allele is the alphabetically last allele at each marker (G, T and G here). For breed composition, format = "matrix" returns the same values with the ids as row names instead of an id column:

geno_mat = vcf_to_dosage("example.vcf", format = "matrix", verbose = FALSE)
geno_mat
##    snp03 snp05 snp09
## M1     0     2     0
## F1     2     0     1
## C1     1     1     1
## attr(,"ploidy")
## [1] 2

See ?vcf_to_dosage and ?ped_to_dosage for the other options, such as reading a dosage field (DS), call-rate filters and multiallelic markers.

3.2 Missing data

There are two kinds of missing information, coded differently:

  • Unknown parent (pedigree): 0. In a fill_pedigree progeny table NA or an empty cell also work.
  • Missing genotype call: NA in a dosage table. The converters turn a VCF call containing . (for example ./.) and a .ped allele coded 0 into NA.

Missing calls never count as errors. Each test only uses the markers where every individual it compares has a call (the tests are explained in How the Mendelian error is calculated):

  • Trio test (validate_pedigree(), and find_parentage() with best_pair or fill_pedigree): markers where the offspring and both parents have a call, the testable markers. A trio with fewer testable markers than min_markers (default 10) is low_markers; one with none is no_data.
  • Single-parent test (recommended corrections, and the single-parent methods of find_parentage()): markers where the parent and the offspring are both homozygous. If there are none, the parent counts as failing (see C_low in Failed trios and recommended correction).
  • Marker summary: mismatch_pct is NA for a marker not tested in any trio.
  • Breed composition: a sample’s missing calls are skipped marker by marker. A marker missing in every reference individual of a line has no frequency (NaN in the output of allele_freq_poly()), which makes solve_composition_poly() stop with an error. Find these markers with rownames(freq)[rowSums(is.na(freq)) > 0] and remove them before rerunning.

4 Pedigree cleaning

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

Issue Meaning Action
Exact duplicates The identical row appears more than once Always removed
Conflicting trios Same id, different parents in different rows Disagreeing parents set to 0 and rows collapsed, or flagged (correct_conflicting_trios = FALSE)
Missing parents A parent is listed but has no row of its own Always added as a row with unknown parents (0)
Inconsistent sex roles Same individual used as sire and as dam Set to 0 as a parent everywhere, or flagged (correct_inconsistent_sex_roles = FALSE)
Cycles / dependencies An individual is its own ancestor, or is recorded as its own parent Flagged only; must be fixed manually, since there is no safe automatic choice. Each cycle is listed once, as a path in which every id is a recorded parent of the one before it (A -> B -> A)

Unknown parents must be coded 0. With verbose = TRUE the function prints a report of every check:

clean_ped_results = check_ped(pedigree, verbose = TRUE)
## 
## === Pedigree Quality Check Report ===
## 
## Exact duplicate trios detected (removed in corrected pedigree):
##    row_number        id male_parent female_parent
## 10         10 C_missing          M1            F1
## 11         11 C_missing          M1            F1
## 
## No conflicting trios found.
## 
## No missing parents found.
## 
## No IDs found as both male_parent and female_parent.
## 
## No dependencies detected.
## 
## The corrected pedigree is included in the returned list as corrected_pedigree.

The result is a list with one table per issue plus the corrected pedigree. An empty table means the issue was not found:

names(clean_ped_results)
## [1] "exact_duplicates"       "conflicting_trios"     
## [3] "inconsistent_sex_roles" "missing_parents"       
## [5] "dependencies"           "corrected_pedigree"
sapply(clean_ped_results, nrow)   # number of rows in each table
##       exact_duplicates      conflicting_trios inconsistent_sex_roles 
##                      2                      0                      0 
##        missing_parents           dependencies     corrected_pedigree 
##                      0                      0                     10
Element In the toy data
exact_duplicates Both copies of the repeated C_missing row; row_number is their position in the input
conflicting_trios Empty: no id has two different sets of parents
inconsistent_sex_roles Empty: no individual is used as both male and female parent
missing_parents Empty: every parent has its own row
dependencies Empty: no cycles
corrected_pedigree The pedigree without the duplicate, ready for validate_pedigree()
clean_ped_results$exact_duplicates
# the corrected, analysis-ready pedigree
clean_ped = clean_ped_results$corrected_pedigree
clean_ped

5 Pedigree validation

validate_pedigree() tests every trio (offspring, male parent, female parent) of the pedigree against the genotypes. It returns the status of each trio, a recommended correction when a trio fails, a corrected pedigree and a per-marker summary. Here the founders are passed with founders_file, and marker_trio_table = TRUE adds the marker-by-trio table used below. Click the label under the code to see the console report.

ped_validate_results = validate_pedigree(clean_ped, genotypes,
                                         founders_file     = founders,
                                         marker_trio_table = TRUE,
                                         verbose = TRUE, plot_results = TRUE)
Show the report printed with verbose = TRUE
## Found 1 trios with missing genotype data; flagged as no_genotype_data.
## 
## --- Trio Validation Summary ---
## Total trios in pedigree: 10
## fail                        : 1 (10.0%)
## founders                    : 4 (40.0%)
## low_markers                 : 1 (10.0%)
## no_genotype_data            : 1 (10.0%)
## pass                        : 3 (30.0%)
## Error threshold: 5%
## Homozygous threshold: 2%
## Minimum markers required: 10
## 
## Correction summary:
##   low_markers_remove_female_parent: 1
##   remove_male_parent: 1
## 
##            id orig_male_parent orig_female_parent
##        <char>           <char>             <char>
##  1:        M1                0                  0
##  2:        F1                0                  0
##  3:        M2                0                  0
##  4:        F2                0                  0
##  5:        C1               M1                 F1
##  6:        C2               M2                 F2
##  7:        C3               M1                 F2
##  8:     C_bad               M2                 F1
##  9:     C_low               M1                 F1
## 10: C_missing               M1                 F1
##     trio_mendelian_error_pct trio_markers_tested           status
##                        <num>               <int>           <char>
##  1:                       NA                   0         founders
##  2:                       NA                   0         founders
##  3:                       NA                   0         founders
##  4:                       NA                   0         founders
##  5:                     0.00                  12             pass
##  6:                     0.00                  12             pass
##  7:                     0.00                  12             pass
##  8:                    41.67                  12             fail
##  9:                     0.00                   4      low_markers
## 10:                       NA                   0 no_genotype_data
##               recommended_correction male_parent_hom_error_pct
##                               <char>                     <num>
##  1:                             none                        NA
##  2:                             none                        NA
##  3:                             none                        NA
##  4:                             none                        NA
##  5:                             none                        NA
##  6:                             none                        NA
##  7:                             none                        NA
##  8:               remove_male_parent                     83.33
##  9: low_markers_remove_female_parent                      0.00
## 10:                             none                        NA
##     female_parent_hom_error_pct
##                           <num>
##  1:                          NA
##  2:                          NA
##  3:                          NA
##  4:                          NA
##  5:                          NA
##  6:                          NA
##  7:                          NA
##  8:                           0
##  9:                          NA
## 10:                          NA

The result is a list. full_results has one row per trio, and pass, fail, low_markers, no_genotype_data, founders and missing_parents hold the same rows split by status. corrected_pedigree, marker_summary and marker_trio_table are described below, and plot is the histogram above.

names(ped_validate_results)
##  [1] "pass"               "fail"               "low_markers"       
##  [4] "no_genotype_data"   "founders"           "missing_parents"   
##  [7] "full_results"       "corrected_pedigree" "marker_summary"    
## [10] "marker_trio_table"  "plot"
sapply(Filter(is.data.frame, ped_validate_results), nrow)   # rows in each table
##               pass               fail        low_markers 
##                  3                  1                  1 
##   no_genotype_data           founders    missing_parents 
##                  1                  4                  0 
##       full_results corrected_pedigree     marker_summary 
##                 10                 10                 12 
##  marker_trio_table 
##                 12
ped_validate_report = ped_validate_results$full_results
ped_validate_report

In full_results, trio_mendelian_error_pct is the trio error over trio_markers_tested testable markers, and male_parent_hom_error_pct and female_parent_hom_error_pct are the single-parent errors, computed only for failed and low-marker trios. The thresholds used here are the defaults; they and the other settings are arguments of the function (see ?validate_pedigree).

5.1 How the Mendelian error is calculated

The trio error uses every marker where the offspring and both parents have a call, heterozygous or homozygous, not only the homozygous ones. At each testable marker the question is whether the two parents can produce the offspring. Each parent passes on a gamete with half of its alleles (ploidy / 2), so a parent with dosage g passes between max(0, g − ploidy / 2) and min(g, ploidy / 2) copies of allele B. The offspring dosage must lie between the sum of the two parents’ minimums and the sum of their maximums; otherwise the marker is a Mendelian error. For a diploid this gives the familiar rules:

Male parent Female parent Possible offspring
0 0 0
0 1 0 or 1
0 2 1
1 1 0, 1 or 2
1 2 1 or 2
2 2 2

For a tetraploid, a parent with dosage 3 passes 1 or 2 copies, so a 3 x 0 cross can only give 1 or 2. The trio error is then

trio error (%) = 100 × markers with an error / testable markers

and the trio passes when it is at or below trio_error_threshold (default 5%). For odd ploidy (e.g. triploid) gametes are unbalanced, so the trio test falls back to the opposite-homozygote check of the single-parent test (see Failed trios and recommended correction): a marker is an error when the offspring is homozygous and a homozygous parent carries the other allele.

Here is the trio of C_bad (recorded as M2 x F1) marker by marker, with the verdict from marker_trio_table:

c_bad_trio = data.frame(
  marker  = ped_validate_results$marker_trio_table$marker,
  M2      = unlist(genotypes[genotypes$id == "M2",    -1]),
  F1      = unlist(genotypes[genotypes$id == "F1",    -1]),
  C_bad   = unlist(genotypes[genotypes$id == "C_bad", -1]),
  verdict = ped_validate_results$marker_trio_table$C_bad_M2_F1,
  row.names = NULL
)
c_bad_trio

5 of the 12 testable markers are errors, a trio error of 41.67%, so the trio fails. At snp01, for example, M2 is 2 and F1 is 0, so every child must be 1, but C_bad is 0. Heterozygous markers count too: at snp10 the cross 2 x 1 can only give 1 or 2, and C_bad is 0. All errors are labelled mismatch_male, meaning M2 cannot be the parent whatever F1 contributes; the single-parent test below confirms this.

5.2 Status of each trio

Status Meaning
pass Trio error at or below trio_error_threshold
fail Trio error above the threshold
low_markers Fewer than min_markers testable markers
no_genotype_data At least one of the three individuals is not in the genotype table
founders Both parents 0 and the id is listed in founders_file
missing_both_parents, missing_male_parent, missing_female_parent One or both parents are 0, and the individual is not a listed founder
no_data No testable marker at all (only in full_results)

In the toy data C1, C2 and C3 pass, C_bad fails, C_low is low_markers and C_missing has no_genotype_data. M1, F1, M2 and F2 are founders because they are listed in founders; without founders_file they would be missing_both_parents.

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

With marker_trio_table = TRUE (as in the run above), the result includes the verdict of every marker for every trio with both parents known and genotyped. Columns are named id_male_female; the C_bad_M2_F1 column was used above.

ped_validate_results$marker_trio_table
Label Meaning
match The trio is consistent at this marker
mismatch_male, mismatch_female Only that parent is inconsistent with the offspring
mismatch_both Each parent alone is inconsistent with the offspring
mismatch_combination Each parent alone is compatible, but the pair cannot produce the offspring
missing_male, missing_female, missing_progeny That individual has no call; every missing individual is named (for example missing_progeny_male_female)
untestable No call is missing, but the marker cannot be tested (odd ploidy, offspring not homozygous)

6 Parentage assignment

When parentage is unknown, find_parentage() tests every candidate in parents against each individual in progeny and reports the best one. How candidates are scored depends on the method; the score is reported as mendelian_error_pct, and lower is better:

Method Candidates tested Score
best_pair (default) Every male x female pair Trio error, as in validate_pedigree(): all markers where the three individuals have a call
best_male_parent, best_female_parent Each candidate of that sex alone Single-parent error: only markers homozygous in both candidate and progeny
best_match Each candidate of either sex alone Single-parent error, as above
fill_pedigree The missing parent(s) of a pedigree; a known parent is kept fixed Trio error (see Filling in missing parents of a pedigree)

Only the single-parent methods rely on homozygous markers. Prefer best_pair when both parents are among the candidates: it uses all markers and both parents together, so it separates candidates much better. Use the single-parent methods when the other parent is not genotyped.

The best candidate is always reported, with a status: pass (score at or below error_threshold, default 5%), high_error (above it: no candidate fits well, so the true parent may be missing from parents), low_markers (fewer than min_markers markers tested) or, for fill_pedigree, known_parent_no_genotype. See ?find_parentage for all options.

find_parentage_results = find_parentage(genotypes, parents, progeny,
                                        method = "best_pair",
                                        verbose = TRUE, plot_results = TRUE)
Show the report printed with verbose = TRUE
## Parent selfing is disallowed. Pairs with identical parents are removed.
## 
## === Parentage Assignment Report ===
## 
## Total progeny evaluated: 5
## Method: best_pair  |  Error threshold: 5%  |  Minimum markers: 10
## 
## 4 progeny passed (80.0%).
## 
## No progeny flagged for high error.
## 
## 1 progeny flagged low_markers (20.0%):
##        id male_parent female_parent mendelian_error_pct markers_tested
##    <char>      <char>        <char>              <char>          <int>
## 1:  C_low          M1            F1                0.00              4
##         status tied_candidates
##         <char>           <int>
## 1: low_markers               1
## 
## Full results are included in the returned list as $full_results.

The result is a list like the one from validate_pedigree(): full_results has one row per progeny, pass, high_error and low_markers hold the same rows split by status, and plot is the histogram above.

names(find_parentage_results)
## [1] "pass"         "high_error"   "low_markers"  "full_results"
## [5] "plot"
sapply(Filter(is.data.frame, find_parentage_results), nrow)   # rows in each table
##         pass   high_error  low_markers full_results 
##            4            0            1            5
find_parentage_report = find_parentage_results$full_results
find_parentage_report

Each row gives the best pair (male_parent, female_parent), its error (mendelian_error_pct), the number of markers tested, the status and tied_candidates, the number of pairs that share the best error. C1, C2 and C3 get their recorded parents at 0% error, and C_bad, which failed validation, is recovered to its true pair (M1 x F1). C_low also gets M1 x F1, but on only 4 markers, so it stays low_markers. C_missing is not in the genotype table, so it is dropped with a warning.

6.1 Top N and ties

tied_candidates above 1 means several candidates share the best error, so the assignment is ambiguous; this is common with few markers or with clones. To see how close the alternatives are, set top_n:

top3_results = find_parentage(genotypes, parents, progeny, method = "best_pair",
                              top_n = 3, verbose = FALSE, plot_results = FALSE)
names(top3_results)
## [1] "pass"         "high_error"   "low_markers"  "full_results"
## [5] "plot"
top3 = top3_results$full_results
top3

The list has the same elements; only full_results gains columns. Rank 1 stays in the base columns (male_parent, female_parent, mendelian_error_pct, markers_tested), and ranks 2 and 3 are added with the suffixes _2 and _3 (male_parent_2, female_parent_2, mendelian_error_pct_2, markers_tested_2, and so on). status and tied_candidates refer to rank 1. A large gap between ranks means a clear assignment: for C_bad the best pair has 0.00% error and the second best (M1 x F2) 25%. With the default top_n = 1, tied best pairs are added in the same suffix columns (show_ties = TRUE).

6.2 Filling in missing parents of a pedigree

With method = "fill_pedigree", progeny is a pedigree (id, male_parent, female_parent), typically the corrected_pedigree from validate_pedigree(). Parents coded 0, NA or empty are unknown. Each row is handled on its own: a known male parent searches only females, a known female parent only males, and a row with both parents unknown gets the full pair search, unless it is a listed founder. Rows with both parents known are skipped. The parents table (with sex) is still required.

filled_results = find_parentage(genotypes, parents,
                                ped_validate_results$corrected_pedigree,
                                method        = "fill_pedigree",
                                founders_file = founders,
                                verbose = FALSE, plot_results = FALSE)
names(filled_results)
## [1] "pass"         "high_error"   "low_markers"  "full_results"
## [5] "plot"
filled_results$full_results

The list has the same elements as before. Only the rows with a parent removed by validation are searched: the founders are skipped because they are in founders, and the complete trios because both parents are known. full_results gains two columns: search_mode (known_male, known_female, or pair when both parents are unknown) and known_parent_error_pct, the single-parent error of the parent that was kept (a high value suggests that it is wrong too). C_bad kept F1, so only males were tested, and M1 is found at 0%. C_low kept M1 and gets F1 back, but on too few markers, so it stays low_markers.

7 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

R2 measures how well the estimated mixture of reference lines explains a sample’s genotypes. The sample’s dosages are divided by the ploidy (0, 0.25, …, 1 for a tetraploid); the predicted value at each marker is the allele frequency of the reference lines weighted by the estimated proportions; and R2 is the squared correlation between the two across the sample’s non-missing markers. It describes the fit, not the accuracy of the proportions. It is usually well below 1 even for a pure individual, because each genotype is one of a few discrete values while the prediction is a frequency, so every marker at which the line is polymorphic adds unexplained variation. Compare R2 between samples analysed with the same reference panel and markers: a sample with a much lower R2 than the rest may come from a population that is not in the reference panel, or have poor-quality genotypes.

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

8 Complete workflow

# genotypes: a dosage table, or pass the VCF / .ped path straight to each function
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, founders_file = founders, ploidy = 4)

# look for the parents removed by validation (and any other missing parent)
filled = find_parentage(genotypes, parents, validated$corrected_pedigree,
                        method = "fill_pedigree", founders_file = founders, ploidy = 4)

# or assign parents to new progeny
assigned = find_parentage(genotypes, parents, progeny, method = "best_pair", 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)

9 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.
  • List your founders in founders_file, so their unknown parents are not reported as missing.
  • 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.

10 How to cite

10.1 BIGpopA

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

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

10.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

11 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.1
## 
## loaded via a namespace (and not attached):
##  [1] tidyselect_1.2.1    viridisLite_0.4.3   dplyr_1.2.1        
##  [4] farver_2.1.2        S7_0.2.2            fastmap_1.2.0      
##  [7] janitor_2.2.1       xopen_1.0.1         rex_1.2.2          
## [10] digest_0.6.39       timechange_0.4.0    lifecycle_1.0.5    
## [13] cluster_2.1.8.3     ellipsis_0.3.3      processx_3.9.0     
## [16] magrittr_2.0.5      compiler_4.5.2      rlang_1.3.0        
## [19] sass_0.4.10         tools_4.5.2         yaml_2.3.12        
## [22] data.table_1.18.6.1 knitr_1.52          labeling_0.4.3     
## [25] prettyunits_1.2.0   pkgbuild_1.4.8      curl_8.0.0         
## [28] RColorBrewer_1.1-3  pkgload_1.5.3       memuse_4.2-3       
## [31] covr_3.6.5          withr_3.0.3         purrr_1.2.2        
## [34] desc_1.4.3          grid_4.5.2          urlchecker_2.0.0   
## [37] ggplot2_4.0.3       scales_1.4.0        MASS_7.3-66        
## [40] cli_3.6.6           rmarkdown_2.32      vegan_2.7-6        
## [43] generics_0.1.4      otel_0.2.0          rstudioapi_0.19.0  
## [46] httr_1.4.9          sessioninfo_1.2.4   ape_5.8-1          
## [49] cachem_1.1.0        stringr_1.6.0       splines_4.5.2      
## [52] parallel_4.5.2      vctrs_0.7.3         devtools_2.5.2     
## [55] Matrix_1.7-6        jsonlite_2.0.0      callr_3.8.0        
## [58] rcmdcheck_1.4.0     checkhelper_1.0.1   testthat_3.3.2     
## [61] vcfR_1.16.0         jquerylib_0.1.4     spelling_2.3.2     
## [64] glue_1.8.1          lubridate_1.9.5     stringi_1.8.9      
## [67] gtable_0.3.6        quadprog_1.5-8      tibble_3.3.1       
## [70] pillar_1.11.1       htmltools_0.5.9     brio_1.1.5         
## [73] R6_2.6.1            pinfsc50_1.3.0      rprojroot_2.1.1    
## [76] evaluate_1.0.5      lattice_0.23-1      memoise_2.0.1      
## [79] snakecase_0.11.1    bslib_0.12.0        Rcpp_1.1.2         
## [82] nlme_3.1-171        permute_0.9-10      mgcv_1.9-4         
## [85] xfun_0.61           fs_2.1.0            usethis_3.2.2      
## [88] pkgconfig_2.0.3