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:
vcf_to_dosage(),
ped_to_dosage()check_ped()validate_pedigree()find_parentage()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):
vcf_to_dosage() and
ped_to_dosage() show the converted table (see Converting genotype files).BIGpopA does not call dosages; use your pipeline first (we recommend BIGapp).
BIGpopA requires R >= 4.4.0. Install it from CRAN by copying and pasting this code into your R console:
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)
pedigreeC_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.
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_ped = ped_to_dosage("example.ped", verbose = FALSE) # example.map in the same folder gives the marker names
geno_pedBoth have the layout of genotypes and, here, the same
values:
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:
## 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.
There are two kinds of missing information, coded differently:
0. In a
fill_pedigree progeny table NA or an empty
cell also work.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):
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.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).mismatch_pct is
NA for a marker not tested in any trio.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.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:
##
## === 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:
## [1] "exact_duplicates" "conflicting_trios"
## [3] "inconsistent_sex_roles" "missing_parents"
## [5] "dependencies" "corrected_pedigree"
## 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() |
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)## 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.
## [1] "pass" "fail" "low_markers"
## [4] "no_genotype_data" "founders" "missing_parents"
## [7] "full_results" "corrected_pedigree" "marker_summary"
## [10] "marker_trio_table" "plot"
## 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
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).
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_trio5 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.
| 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.
A failed trio shows that the pair is wrong, but not which parent. Each parent is therefore tested alone with the single-parent test, which uses only the markers where the parent and the offspring are both homozygous. A parent homozygous for one allele passes that allele to every child, so a child homozygous for the other allele (0 vs 2 in a diploid) cannot be its offspring. Heterozygous markers are left out because, without the other parent, they rarely rule a parent out.
single-parent error (%) = 100 × markers where parent and offspring are homozygous for different alleles / markers where both are homozygous
A parent fails when this error is above
single_parent_error_threshold (default 2%), or when there
is no marker to compare. For C_bad, M2 has a
single-parent error of 83.33% and F1 of 0%, so the
recommendation is remove_male_parent. The
recommended_correction column gives the result for every
trio:
recommended_correction |
Interpretation |
|---|---|
none |
Trio passes, or it is a founder, has a missing parent, has no genotype data or has no testable marker |
remove_male_parent |
The trio fails and only the male parent fails the single-parent test |
remove_female_parent |
The trio fails and only the female parent fails the single-parent test |
remove_both |
The trio fails and both parents fail the single-parent test |
unresolved |
The trio fails but neither parent fails the single-parent test, for example an offspring recorded as its own parent; the pedigree is probably wrong but the faulty parent cannot be identified |
low_markers_* |
The trio has fewer than min_markers testable markers;
the suffix is the decision the same rules would give
(low_markers_keep_both means the error is within the
threshold and both parents pass) |
C_low shows the weak point of having few markers. Its
trio error is 0%, but F1 is heterozygous at all 4 of its
testable markers, so the single-parent test has nothing to compare for
F1. F1 then counts as failing, and the
recommendation is low_markers_remove_female_parent. Treat
low_markers_* recommendations with caution, and genotype
these samples again if possible.
corrected_pedigree sets the parent(s) named by the
recommendation to 0, here the male parent of
C_bad and the female parent of C_low.
unresolved trios are left unchanged for manual review. Run
find_parentage() with method = "fill_pedigree"
on the corrected pedigree to look for the right parents (see Filling in missing
parents of a pedigree).
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.
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.
| 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) |
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)## 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.
## [1] "pass" "high_error" "low_markers" "full_results"
## [5] "plot"
## pass high_error low_markers full_results
## 4 0 1 5
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.
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"
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).
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"
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.
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))
predictionR2 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_resultsSample1 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).
# 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)check_ped() ->
validate_pedigree() -> find_parentage()
(fill_pedigree) on the corrected pedigree.founders_file, so their unknown
parents are not reported as missing.unresolved flag.tied_candidates.marker_summary before blaming individuals; remove
bad markers and rerun.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.
## 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