Genotype quality control

Filters animals and markers by call rate, markers by minor allele frequency and (optionally) a Hardy-Weinberg chi-square, and, with pedigree =, the offspring whose genotype contradicts a genotyped parent; and REPORTS what each filter removed: silent filtering is how a panel quietly loses the markers that mattered.

Usage

qc_genotypes(
  m,
  min_call_rate = 0.9,
  min_maf = 0.01,
  hwe_p = NULL,
  min_call_rate_animal = getOption("BreedingR.min_call_rate_animal", 0.9),
  pedigree = NULL,
  exclusion_threshold = 0.01,
  sex_markers = if (is.list(m)) m[["sex_markers", exact = TRUE]]
)

Arguments

m

dosage matrix 0/1/2/NA (animals in rows), or a genotypes = list(ids, m); with the list the result carries ids too and goes straight to genotypes =. Double, integer and raw storage give the same filtering and counts; in a raw matrix, as read_blupf90_snp() returns it, 5 is the missing code

min_call_rate

keep markers with at least this fraction of calls among the kept animals

min_maf

keep markers with minor allele frequency at least this

hwe_p

if not NULL, drop markers whose HWE chi-square p-value falls below it

min_call_rate_animal

keep animals with at least this fraction of calls over the kept markers that vary; the default is the line of the single-step fitters, getOption("BreedingR.min_call_rate_animal", 0.90). NULL keeps every animal, and then only the markers are filtered

pedigree

NULL, or a pedigree (animal, sire, dam): then the offspring in conflict with a genotyped parent are dropped. Needs m as a list(ids, m)

exclusion_threshold

fraction of opposite homozygotes, among the markers called in both, above which a parent-offspring pair is a conflict

sex_markers

the columns of m on the sex chromosomes (indices, names or a logical vector), left out of the parentage test only; the map of read_plink() tells them apart. By default m$sex_markers when m is a list, so a second qc_genotypes() on a result reads what the first one was given

Value

list with m (filtered matrix, in the storage it came in: a raw matrix stays raw, with 5 as missing), ids (the kept ids, from the list or the rownames; NULL without either), rows and columns (the kept rows and columns of the input), the counts n_animal (animals removed by call rate), n_call, n_maf, n_hwe (markers removed by each filter) and n_parentage (offspring dropped), passes, removed (one row per removed animal: row, id, reason, call_rate over the markers kept at its pass, and pass), with pedigree =, parentage (one row per genotyped pair, as parentage_check() gives it, from the last pass that tested it) and, with sex_markers, sex_markers (one logical per column of the returned m). Every fitter reads genotypes$sex_markers and leaves those columns out of its parentage warning, so what passed the parentage filter here is not flagged again in the fit

Details

The call rates are filtered worst first, in passes. Each pass counts again the call rate of every kept animal (its calls over the kept markers that vary, the ones that enter G) and of every kept marker (its calls over the kept animals); the side holding the lowest rate loses its entries below the threshold that are also below the lowest failing rate of the other side. A failed sample thus leaves without taking the markers with it, and a bad batch of markers leaves without taking the animals. PLINK applies --mind before --geno in a single pass, which drops every animal when a batch of markers failed; the usual protocols filter markers first and run a lenient pass before the strict one for the same reason (Marees et al. 2018), and the passes here are that idea carried to the end. With two panels the side that misses more goes: the animals of a panel holding a small part of the markers (impute them first to keep them), or the markers that most animals lack (the shared ones stay). Then MAF and HWE on the kept animals, and the parentage test; anything they remove changes the rates, and the passes repeat until nothing changes. So the result always passes the call-rate check of the single-step fitters (see model()), whose line is the same options(BreedingR.min_call_rate_animal) that gives the default here.

Heterogeneous missingness takes many passes (hundreds when the failing animals and markers interleave), but the matrix is read once: the counts are kept and updated as rows and columns leave, reading only the rows and columns that left. A pass then costs a few vector operations over the animals and markers, and the parentage test, with pedigree =, reads the genotyped pairs once each time it runs.

The 0.90 for animals is the callrateAnim default of preGSf90, counted over the markers that vary where preGSf90 counts every SNP, and 0.90 for markers is its callrate. The parentage filter is the verify_parentage 3 default of preGSf90: a genotyped parent and offspring with opposite homozygotes on more than exclusion_threshold of the markers called in both cannot be parent and offspring, and the OFFSPRING's genotype is dropped (preGSf90's exclusion_threshold is the same 1%, over every SNP); parentage_check() lists the pairs, for the alternative of setting the parent unknown in the pedigree instead. Leave the sex-chromosome markers out of that test with sex_markers = (see parentage_check() for why).

A removed animal leaves the genotypes only; it stays in the pedigree. The removals are a message, and a warning when more than half of the animals or of the markers went or when rows leave a bare matrix without rownames (the ids live elsewhere: subset them with ids[result$rows]). Leaving no animal or no marker is an error that says what removed them.

References

Chang, C.C., Chow, C.C., Tellier, L.C.A.M., Vattikuti, S., Purcell, S.M. & Lee, J.J. (2015). Second-generation PLINK: rising to the challenge of larger and richer datasets. GigaScience 4:7 (PLINK 1.9, whose order of operations applies --mind before --geno).

Marees, A.T., de Kluiver, H., Stringer, S., Vorspan, F., Curis, E., Marie-Claire, C. & Derks, E.M. (2018). A tutorial on conducting genome-wide association studies: quality control and statistical analysis. International Journal of Methods in Psychiatric Research 27:e1608 (markers before individuals, a relaxed threshold before the strict).

Misztal, I. et al. (2018). Manual for BLUPF90 family of programs, preGSf90 quality control options (callrate, callrateAnim, verify_parentage, exclusion_threshold, sex_chr).

Examples

s <- simulate_breeding(20, 2, 20, n_markers = 300, seed = 1)
g <- s$genotypes
g$m[3, 1:200] <- NA           # a failed sample
g$m[, 1:30][sample(length(g$m[, 1:30]), 1000)] <- NA   # a bad batch of markers
q <- qc_genotypes(g, pedigree = s$pedigree)
#> qc_genotypes: kept 59 of 60 animal(s) and 269 of 300 marker(s) in 6 pass(es). Animals left out of the genotypes (they stay in the pedigree): 1 with a call rate below 0.90 ('F003': 0.335). Markers left out: 30 by call rate below 0.90, 1 by MAF. The kept rows are in $rows.
q[c("n_animal", "n_call", "n_maf", "passes")]
#> $n_animal
#> [1] 1
#> 
#> $n_call
#> [1] 30
#> 
#> $n_maf
#> [1] 1
#> 
#> $passes
#> [1] 6
#> 
q$removed
#>   row   id    reason call_rate pass
#> 1   3 F003 call rate 0.3355705    3