Package {BIGpopA}

Compute Allele Frequencies for Populations

Description

Computes allele frequencies for specified populations given SNP array data.

Usage

allele_freq_poly(geno, populations, ploidy = 2)

Arguments

Details

When geno is a .ped file, the returned matrix carries a counted_allele attribute. solve_composition_poly() uses it to code a validation .ped file with the same counted allele at each marker.

Value

A matrix of allele frequencies with SNPs in rows and populations in columns.

References

Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP. Estimation of genome-wide and locus-specific breed composition in pigs. Transl Anim Sci. 2017 Feb 1;1(1):36-44.

Examples

geno_matrix <- matrix(
  c(4, 1, 4, 0,
    2, 2, 1, 3,
    0, 4, 0, 4,
    3, 3, 2, 2,
    1, 4, 2, 3),
  nrow = 4, ncol = 5, byrow = FALSE,
  dimnames = list(paste0("Ind", 1:4), paste0("S", 1:5))
)

pop_list <- list(
  PopA = c("Ind1", "Ind2"),
  PopB = c("Ind3", "Ind4")
)

allele_freqs <- allele_freq_poly(geno = geno_matrix,
                                 populations = pop_list,
                                 ploidy = 4)
print(allele_freqs)

Assign Each Sample to a Reference Line

Description

Builds the composition table shown in the Familia BreedTools tab from the output of solve_composition_poly(): one row per sample with the estimated percentage of each reference population and the predicted line. A sample is assigned to the population with its highest proportion only if that proportion is at least assign_threshold; otherwise it is labelled "Undetermined".

Usage

assign_composition_line(composition, assign_threshold = 50)

Arguments

Value

A data.frame with an ID column, one column per reference population holding the estimated percentage (0 to 100, rounded to one decimal) and a ⁠Predicted line⁠ column. Ties between populations go to the first column.

Examples

comp <- rbind(
  Sample01 = c(Duroc = 0.998, Hampshire = 0.000, Landrace = 0.002, R2 = 0.991),
  Sample02 = c(Duroc = 0.529, Hampshire = 0.471, Landrace = 0.000, R2 = 0.973),
  Sample03 = c(Duroc = 0.340, Hampshire = 0.330, Landrace = 0.330, R2 = 0.952)
)
assign_composition_line(comp)
assign_composition_line(comp, assign_threshold = 0)

Check and Correct Common Pedigree Errors

Description

Reads a 3-column pedigree file (id, male_parent, female_parent) and performs quality checks, optionally correcting detected errors. Exact duplicates and missing parents are always corrected. Conflicting trios and inconsistent sex roles are corrected when their respective arguments are TRUE. Cycles are reported only and must be resolved manually.

Usage

check_ped(
  ped.file,
  seed = NULL,
  verbose = TRUE,
  correct_conflicting_trios = TRUE,
  correct_inconsistent_sex_roles = TRUE
)

Arguments

Value

An invisible named list of data frames:

exact_duplicates

Exact duplicate rows found in the input.

conflicting_trios

IDs with conflicting male_parent or female_parent assignments.

inconsistent_sex_roles

Rows where a conflicting ID appears as male_parent or female_parent.

missing_parents

Parent IDs absent from id, added as founders.

dependencies

Cycles detected in the pedigree. Must be resolved manually.

corrected_pedigree

Corrected pedigree table.

Author(s)

Josue Chinchilla-Vargas

Examples

# Self-contained example using a data.frame
ped_df <- data.frame(
  id            = c("A", "B", "C", "C", "D"),
  male_parent   = c("0", "0", "A", "A", "B"),
  female_parent = c("0", "0", "B", "B", "C"),
  stringsAsFactors = FALSE
)
ped_errors <- check_ped(ped.file = ped_df, seed = 101919, verbose = FALSE)
names(ped_errors)
head(ped_errors$corrected_pedigree)


library(data.table)
ped_dt <- data.table(id = c("A", "B", "C"),
                     male_parent   = c("0", "0", "A"),
                     female_parent = c("0", "0", "B"))
ped_errors <- check_ped(ped.file = ped_dt, verbose = FALSE)

Find Parentage Assignments for Progeny

Description

Assigns the most likely parent(s) to each progeny from SNP genotype data using Mendelian error rates or homozygous mismatch rates. Parents or progeny absent from the genotype file are removed with a warning.

Usage

find_parentage(
  genotypes_file,
  parents_file,
  progeny_file,
  method = "best_pair",
  min_markers = 10,
  error_threshold = 5,
  show_ties = TRUE,
  allow_parent_selfing = FALSE,
  exclude_self_match = TRUE,
  verbose = TRUE,
  plot_results = TRUE,
  ploidy = 2,
  top_n = 1,
  founders_file = NULL
)

Arguments

Details

Method "fill_pedigree" ("fill in the blanks"). The progeny_file must have male_parent and female_parent columns. A parent coded 0, NA or empty is treated as unknown and each row is handled on its own: a known male searches only females, a known female searches only males, and a row with both unknown gets the full pair search (unless listed in founders_file). Rows with both parents known are skipped. The search is the "best_pair" Mendelian trio test with the known parent held fixed, and candidates come from parents_file (sex M/F/A as before); the known parent only needs genotypes. Two columns are added to the results: search_mode ("pair", "known_male" or "known_female") and known_parent_error_pct (homozygous mismatch between the known parent and the progeny; a high value suggests the known parent is wrong). A known parent absent from the genotype file gives status = "known_parent_no_genotype". Without method = "fill_pedigree", any male_parent / female_parent columns in progeny_file are ignored.

Value

A named list (returned invisibly) with elements:

pass

Progeny with a confident parentage assignment.

high_error

Progeny whose best assignment exceeds the error threshold.

low_markers

Progeny with insufficient markers for a valid assignment.

full_results

Complete data.table with all progeny and all output columns. tied_candidates gives the number of candidates (pairs for "best_pair") that share the lowest error for each progeny; values above 1 mean the best assignment is ambiguous.

plot

ggplot object if plot_results = TRUE, otherwise NULL.

Author(s)

Josue Chinchilla-Vargas

Examples

geno_df <- data.frame(
  id  = c("P1", "P2", "P3", "Off1", "Off2"),
  S1  = c(0L, 2L, 0L, 1L, 0L),
  S2  = c(2L, 0L, 2L, 1L, 2L),
  S3  = c(0L, 2L, 0L, 1L, 0L),
  S4  = c(2L, 0L, 2L, 1L, 2L),
  S5  = c(0L, 2L, 0L, 1L, 0L),
  S6  = c(2L, 0L, 2L, 1L, 2L),
  S7  = c(0L, 2L, 0L, 1L, 0L),
  S8  = c(2L, 0L, 2L, 1L, 2L),
  S9  = c(0L, 2L, 0L, 1L, 0L),
  S10 = c(2L, 0L, 2L, 1L, 2L)
)

parents_df <- data.frame(
  id  = c("P1", "P2", "P3"),
  sex = c("M",  "F",  "F"),
  stringsAsFactors = FALSE
)

progeny_df <- data.frame(
  id = c("Off1", "Off2"),
  stringsAsFactors = FALSE
)

results <- find_parentage(
  genotypes_file = geno_df,
  parents_file   = parents_df,
  progeny_file   = progeny_df,
  method         = "best_pair",
  verbose        = FALSE,
  plot_results   = FALSE
)
print(results$full_results)

Convert a PLINK .ped File to a BIGpopA Genotype Matrix

Description

Reads a PLINK 1 text pedigree file (.ped) and returns allele dosages in the layout used by the rest of BIGpopA: individuals in rows, markers in columns, coded as the number of copies of the counted allele (0, 1, 2) with NA for missing calls. PLINK .ped files are diploid only.

Usage

ped_to_dosage(
  ped,
  map = "auto",
  format = c("data.frame", "matrix"),
  counted_allele = NULL,
  missing_codes = c("0", "N", "-", "."),
  verbose = TRUE
)

Arguments

Details

The .ped format stores allele letters and has no reference/alternate pair, so the counted allele (allele B) is chosen per marker. By default it is the alphabetically / numerically last allele observed at the marker, so 1/2 coded files count allele 2 and an A/G marker counts G. To code a second file the same way as a first one (for example validation genotypes against a reference panel), pass the first file's counted_allele attribute through the counted_allele argument. Otherwise a marker that is monomorphic in one file can be coded in the opposite direction.

Only the individual ID (column 2) is used from the six leading pedigree columns. A .map file is optional and only supplies marker names; it is not used for allele coding. Without it, markers are named SNP1, SNP2, ... in file order.

Value

A data.table with an id column followed by marker columns (format = "data.frame"), or an integer matrix of individuals by markers (format = "matrix"). The attribute "ploidy" is 2 and the attribute "counted_allele" holds the allele counted at each returned marker. Markers with more than two alleles are dropped with a warning.

Author(s)

Josue Chinchilla-Vargas

See Also

vcf_to_dosage(), find_parentage(), allele_freq_poly()

Examples

ped_path <- tempfile(fileext = ".ped")
writeLines(c(
  "F1 P1 0 0 1 -9 A A G G A G",
  "F1 P2 0 0 2 -9 A G G G G G",
  "F1 O1 P1 P2 1 -9 A G G G A G"
), ped_path)

genotypes <- ped_to_dosage(ped_path, map = NULL, verbose = FALSE)
print(genotypes)
attr(genotypes, "counted_allele")

unlink(ped_path)

Compute Genome-Wide Breed Composition

Description

Computes genome-wide breed/ancestry composition using quadratic programming on a batch of animals.

Usage

solve_composition_poly(Y, X, ploidy = 2)

Arguments

Value

A matrix with one row per animal, one column per reference population (estimated proportions, summing to 1), and an R2 column.

References

Funkhouser SA, Bates RO, Ernst CW, Newcom D, Steibel JP. Estimation of genome-wide and locus-specific breed composition in pigs. Transl Anim Sci. 2017 Feb 1;1(1):36-44.

Examples

allele_freqs_matrix <- matrix(
  c(0.625, 0.500,
    0.500, 0.500,
    0.500, 0.500,
    0.750, 0.500,
    0.625, 0.625),
  nrow = 5, ncol = 2, byrow = TRUE,
  dimnames = list(paste0("SNP", 1:5), c("VarA", "VarB"))
)

val_geno_matrix <- matrix(
  c(2, 1, 2, 3, 4,
    3, 4, 2, 3, 0),
  nrow = 2, ncol = 5, byrow = TRUE,
  dimnames = list(paste0("Test", 1:2), paste0("SNP", 1:5))
)

composition <- solve_composition_poly(Y = val_geno_matrix,
                                      X = allele_freqs_matrix,
                                      ploidy = 4)
print(composition)

Validate Pedigree Trios Using Mendelian Error Analysis

Description

Validates parent-offspring trios against SNP genotype data using Mendelian error rates. Identifies incorrect parentage assignments and outputs a corrected pedigree. Founder trios (both parents coded as 0) are preserved unchanged if a founders file is supplied. Trios absent from the genotype file are retained as no_genotype_data.

Usage

validate_pedigree(
  pedigree_file,
  genotypes_file,
  founders_file = NULL,
  trio_error_threshold = 5,
  min_markers = 10,
  single_parent_error_threshold = 2,
  verbose = TRUE,
  plot_results = TRUE,
  ploidy = 2,
  marker_trio_table = FALSE
)

Arguments

Details

The cells of marker_trio_table take one of these values: match; mismatch_male or mismatch_female (the offspring is impossible with that parent alone, whatever the other parent contributes); mismatch_both (the offspring is impossible with either parent alone); mismatch_combination (each parent is compatible alone but the pair is not, for example a 0 x 0 cross with a heterozygous offspring); untestable (odd ploidy, marker not informative); ⁠missing_*⁠ (no call, named after every individual without a call: missing_progeny, missing_male, missing_female, missing_male_female, missing_progeny_male, missing_progeny_female or missing_progeny_male_female).

The recommended_correction column of full_results takes these values:

none

Trio passes, or a parent is missing / no genotype data.

remove_male_parent, remove_female_parent

The trio fails and only one parent fails the single-parent homozygous check; that parent is set to 0 in corrected_pedigree.

remove_both

The trio fails and both parents fail the homozygous check; both are set to 0.

unresolved

The trio fails but neither parent fails the homozygous check (for example a recorded self whose offspring carries alleles the parent lacks). The pedigree is likely wrong but the faulty parent cannot be identified, so corrected_pedigree is left unchanged for manual review.

keep_both

Only with low markers: the trio is within the threshold and both parents pass.

Low-marker trios get the same values prefixed with low_markers_. In the plot, every trio above the error threshold is coloured as a failure; unresolved trios are shown in grey.

Value

An invisible named list with the following elements:

pass

Trios that passed the Mendelian error threshold.

fail

Trios that failed the Mendelian error threshold.

low_markers

Trios with insufficient markers for evaluation.

no_genotype_data

Trios absent from the genotype file.

founders

Trios identified as founders.

missing_parents

Trios with one or both parents coded as 0 (non-founders).

full_results

Complete data.table with all trios and all output columns.

corrected_pedigree

Pedigree table after applying recommended corrections.

marker_summary

Per-marker Mendelian mismatch counts over every trio with both parents known and genotyped: marker, trios_tested, trios_mismatch and mismatch_pct, sorted worst first. Markers with a much higher rate than the rest (e.g. paralogous or miscalled markers) are candidates for removal before rerunning.

marker_trio_table

If marker_trio_table = TRUE, a data.table with one row per marker (marker) and one column per trio with both parents known and genotyped, named id_male_parent_female_parent; otherwise NULL. See Details for the status values.

plot

ggplot object if plot_results = TRUE, otherwise NULL.

Author(s)

Josue Chinchilla-Vargas

Examples

geno_df <- data.frame(
  id  = c("P1", "P2", "P3", "Off1", "Off2"),
  S1  = c(0L, 2L, 0L, 1L, 0L),
  S2  = c(2L, 0L, 2L, 1L, 2L),
  S3  = c(0L, 2L, 0L, 1L, 0L),
  S4  = c(2L, 0L, 2L, 1L, 2L),
  S5  = c(0L, 2L, 0L, 1L, 0L),
  S6  = c(2L, 0L, 2L, 1L, 2L),
  S7  = c(0L, 2L, 0L, 1L, 0L),
  S8  = c(2L, 0L, 2L, 1L, 2L),
  S9  = c(0L, 2L, 0L, 1L, 0L),
  S10 = c(2L, 0L, 2L, 1L, 2L)
)

ped_df <- data.frame(
  id            = c("Off1", "Off2"),
  male_parent   = c("P1",   "P1"),
  female_parent = c("P2",   "P3"),
  stringsAsFactors = FALSE
)

results <- validate_pedigree(
  pedigree_file  = ped_df,
  genotypes_file = geno_df,
  verbose        = FALSE,
  plot_results   = FALSE
)
print(results$full_results)

Convert a VCF File to a BIGpopA Genotype Matrix

Description

Reads a variant call format (VCF) file, or an existing vcfR object, and returns allele dosages in the layout used by the rest of BIGpopA: individuals in rows, markers in columns, coded as the dosage of the alternate (B) allele (0, 1, ..., ploidy) with NA for missing calls.

Usage

vcf_to_dosage(
  vcf,
  ploidy = NULL,
  format = c("data.frame", "matrix"),
  dosage_field = NULL,
  ref_dosage = FALSE,
  biallelic_only = TRUE,
  marker_names = c("auto", "id", "position"),
  min_marker_call_rate = 0,
  min_sample_call_rate = 0,
  verbose = TRUE,
  ...
)

Arguments

Details

The function is ploidy agnostic. Dosage is obtained by counting the non-reference alleles in each GT call, so a diploid 0/1 returns 1 and a tetraploid 0/1/1/1 returns 3, with no assumption about the ploidy level. When ploidy is not supplied it is inferred from the number of allele slots in the calls themselves (the most common slot count across the file). Calls whose slot count disagrees with the file ploidy, for instance a diploid call in an otherwise tetraploid file, are set to missing and reported, since the downstream functions require a single ploidy level.

Phased (|) and unphased (/) separators are both accepted. A call is treated as missing when any of its allele slots is ..

vcfR is used to read the file and is listed under Suggests, so it must be installed before calling this function (install.packages("vcfR")).

Value

A data.table with an id column followed by marker columns (format = "data.frame"), or a numeric matrix of individuals by markers (format = "matrix"). The attribute "ploidy" records the ploidy used.

Author(s)

Josue Chinchilla-Vargas

See Also

find_parentage(), validate_pedigree(), allele_freq_poly()

Examples

if (requireNamespace("vcfR", quietly = TRUE)) {
  # Write a small tetraploid VCF (3 markers x 4 samples) to a temporary file
  vcf_path <- tempfile(fileext = ".vcf")
  writeLines(c(
    "##fileformat=VCFv4.3",
    "##FORMAT=<ID=GT,Number=1,Type=String,Description=\"Genotype\">",
    paste("#CHROM", "POS", "ID", "REF", "ALT", "QUAL", "FILTER", "INFO",
          "FORMAT", "S1", "S2", "S3", "S4", sep = "\t"),
    paste("chr1", "100", "snp1", "A", "T", ".", ".", ".", "GT",
          "0/0/0/0", "0/0/0/1", "0/1/1/1", "1/1/1/1", sep = "\t"),
    paste("chr1", "200", "snp2", "C", "G", ".", ".", ".", "GT",
          "0/0/1/1", "./././.", "0/0/0/0", "0/1/1/1", sep = "\t"),
    paste("chr1", "300", "snp3", "G", "A", ".", ".", ".", "GT",
          "1/1/1/1", "1/1/1/1", "0/0/1/1", "0/0/0/1", sep = "\t")
  ), vcf_path)

  # tetraploid VCF -> id + marker columns for find_parentage()
  genotypes <- vcf_to_dosage(vcf_path, ploidy = 4)
  print(genotypes)

  # same data as an individuals x markers matrix for allele_freq_poly()
  geno_mat <- vcf_to_dosage(vcf_path, ploidy = 4, format = "matrix")

  unlink(vcf_path)
}