Skip to contents

Filter markers based on Hardy-Weinberg Equilibrium

Source: R/filter_hwe.R
filter_hwe.Rd

Testing markers for Hardy-Weinberg proportions is a valuable tool for the analysis and quality control of RADseq datasets. HWE can highligh genotyping errors, presence of null alleles, sequence duplication, copy number variation and other sequencing problems related to read depth. This function is designed for bi-allelic markers and uses the exact test from the package HardyWeinberg. The function is speedy because it uses the C++ code developed by Christopher Chang and also available in PLINK. The p-value generated by the function is the mid p-value. It's computed as half the probability of the current sample + the probabilities of all samples that are more extreme (see references below). Several output are generated to help users filter the data (see details).

sampling sites vs well defined populations: be careful what strata you're investigating and adjust filtering threshold accordinghly for downstream analysis. If you're still at the populations discovery steps, markers under HWD are totally normal.

Prior to HW filtering, I highly recommend removing outlier individuals, filtering coverage and genotype likelihood (see details).

Filter target: Markers.

Usage

filter_hwe(
  data,
  interactive.filter = TRUE,
  filter.hwe = TRUE,
  strata = NULL,
  hw.pop.threshold = NULL,
  midp.threshold = 4L,
  filename = NULL,
  parallel.core = parallel::detectCores() - 1,
  verbose = TRUE,
  ...
)

Arguments

data

A GDS file or connection, or a supported tidy genomic object. Use read_genome to import genomic files.

interactive.filter

Logical indicating whether an interactive filtering session may display diagnostics and ask for thresholds. Default: interactive.filter = TRUE.

filter.hwe

(optional, logical) Used inside radr pipeline. Default: filter.hwe = TRUE.

strata

Optional strata file or object containing individual and group information. Default: strata = NULL.

hw.pop.threshold

(integer, optional) Remove markers that have a certain number of pops in Hardy-Weinberg disequilibrium. With default, all populations in dataset need to be in HWD before discarding the marker. Default: hw.pop.threshold = NULL.

midp.threshold

(integer, optional) By default the function generates blacklists/whitelists of markers and filtered tidy datasets for the 5 mid p-value. However, to get a final filtered object associated with the output of the function, user need to choose one of the 5 mid p-value:

  • 1 = 0.05

  • 2 = 0.01

  • 3 = 0.001

  • 4 = 0.0001

  • 5 = 0.00001

. With default, a very conservative mid p-value threshold is selected. Default: midp.threshold = 4L.

filename

Optional output prefix reserved for compatibility. HWE candidate datasets currently use standardized descriptive filenames in the function output folder. Default: filename = NULL.

parallel.core

Number of workers available for parallel operations. Default: parallel.core = parallel::detectCores() - 1.

verbose

Logical indicating whether progress messages are emitted. Default: verbose = TRUE.

...

Additional arguments passed to lower-level screening or filtering functions.

Value

The filtered data. For GDS workflows, marker filtering metadata is synchronized with the selected HWE result. HWE summaries, plots, candidate datasets, marker lists, and filtering parameters are written to the output folder so threshold sensitivity can be reviewed.

Written in the folder:

  1. genotypes.summary.tsv: A tibble with these columns: MARKERS, STRATA, HET, HOM_ALT, HOM_REF, MISSING, N, FREQ_ALT, FREQ_REF, FREQ_HET, FREQ_HOM_REF_O, FREQ_HET_O, FREQ_HOM_ALT_O, FREQ_HOM_REF_E, FREQ_HET_E, FREQ_HOM_ALT_E, N_HOM_REF_EXP, N_HET_EXP, N_HOM_ALT_EXP, HOM_REF_Z_SCORE, HOM_HET_Z_SCORE, HOM_ALT_Z_SCORE, READ_DEPTH

  2. hw.pop.sum.tsv: a summary tibble with populations, number of markers in total, number of markers monomorphic for the populations, number of markers in Hardy-Weinberg Equilibrium (HWE), number of markers in Hardy-Weinberg Dquilibrium (HWD) with all the different mid p-values observed on the data.

  3. hwd.helper.table.tsv: useful tibble that highlight the number of markers blacklisted based on the number of populations in HWD and mid p-value thresholds.

  4. hwd.plot.blacklist.markers.pdf: useful figure that highlight the number of markers blacklisted based on the number of populations in HWD and mid p-value thresholds.

  5. hwe.manhattan.plot.pdf: manhattan plot of markers in Hardy-Weinberg disequilibrium.

  6. tidy.filtered.hwe.xxx.mid.p.value.xxx.hw.pop.threshold.arrow.parquet: candidate datasets generated for different mid-p and population thresholds

  7. whitelist.markers.hwe.xxx.mid.p.value.xxx.hw.pop.threshold.tsv: several whitelist of markers with different mid p value and populations in HWD thresholds

  8. blacklist.markers.hwd.xxx.mid.p.value.xxx.hw.pop.threshold.tsv: several blacklist of markers with different mid p value and populations in HWD thresholds

Details

HWE threshold

I recommend starting with a low threshold. Serious genotyping errors will generate extreme p-values (e.g. 1e-50), which are detected by any reasonable configuration of this test, while various life-history caracteristics will deflate/inflate Hardy-Weinberg equilibrium. Consequently, it's dangerous to choose a threshold that filters out too many markers.

strategies: Disk space is cheap! Consequently, the function will automatically generate several blacklists/whitelists of markers and filtered tidy data (in the directory) based on the hw.pop.threshold for 5 groups of mid p-values:

  • 1: MID_P_VALUE <= 0.00001

  • 2: MID_P_VALUE <= 0.0001

  • 3: MID_P_VALUE <= 0.001

  • 4: MID_P_VALUE <= 0.01

  • 5: MID_P_VALUE <= 0.05

Test the sensitivity of downstream analysis and delete unwanted datasets.

missing data

The mid-p adjustment tends to bring the null rejection rate in line with the nominal p-value, and also reduces the filter's tendency to favor retention of SNPs with missing data (Graffelman and Moreno, 2013, Purcell et al., 2007).

If pattern of missing data is present in the dataset, or when missing data accross markers vary by more than 0.10, you should not apply a single mid-p-value threshold accross markers.

Read depth, pooling lanes/chips and weird pattern of individual heterozygosity

Because of read depth, heterozygote deficiency is usually observed in RADseq data, but if sequencing lanes/chips were combined to generate individuals with more coverage, the situation will likely be the reverse: heterozygote excess. If lanes/chips pooling was used or if highly variable sequencing coverage is observed between individuals and/or markers, there's a couple qc and filtering steps you should do before conducting HWE filtering.

  • run radr detect_mixed_genomes and detect_het_outliers to highlight outlier individuals with potential heterozygosity problems and get an idea of the genotyping and heterozygote miscall rate

  • I recommend normalizing the data before de novo assembly, if this is not possible...

  • use radr::explore_genomes or a more chirurgical approach to coverage and genotype likelihood with radr::explore_genomes.

Permutation test Hardy-Weinberg equilibrium refers to the statistical independence of alleles within individuals. This independence can also be assessed by permutation test inside the HardyWeinberg package of Jan Graffelman. To filter out markers with genotyping problems the approach provided in this function is enough.

Power test for HWE Based on allele count/frequency and sample size. This function is longer to generate for each markers and is on my todo list to include it in this filter.

Markers under selection and genome scans:

Scared of deleting those precious markers or that the filter might interfere with genome scan analysis/detection ? Don't be. Your markers or analysis is no good if it's done on bad data... Test the sensitivity of your downstream analysis with the datasets generated with the different thresholds.

Note

Hardy-Weinberg assumptions (refresh):

  1. Diploid organisms

  2. Only reproduction sexual occurs

  3. Generations are non-overlapping

  4. Mating is random

  5. Population size is infinitely large

  6. Allele frequencies are equal in the sexes

  7. Migration, Mutation and Selection are negligible

Interactive version

The function calculates HWE results, writes the diagnostic tables and plots, and asks the user to inspect them before filtering. The questions are:

  1. "Do you want to continue with the filtering? (y/n):"

  2. If yes, "Based on figures and tables enter the hw.pop.threshold (integer):"

  3. Finally, select one mid-p threshold: 1 = 0.05, 2 = 0.01, 3 = 0.001, 4 = 0.0001, or 5 = 0.00001.

The population threshold controls how many strata may show disequilibrium before a marker is blacklisted. Answering no retains the diagnostics without applying an HWE filter. Use interactive.filter = FALSE with explicit hw.pop.threshold and midp.threshold values for reproducibility.

References

Weir, B.S. (1996) Genetic data analysis II. Sinauer Associates, Massachusetts. See Chapter3.

Wigginton, J.E., Cutler, D.J. and Abecasis, G.R. (2005) A note on exact tests of Hardy-Weinberg equilibrium, American Journal of Human Genetics (76) pp. 887-893.

Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M.A.R., Bender, D., Maller, J., Sklar, P., de Bakker, P.I.W., Daly, M.J. and Sham, P.C. (2007) PLINK: A Tool Set for Whole-Genome Association and Population-Based Linkage Analyses. American Journal of Human Genetics 81(3) pp. 559-575.

Graffelman, J. and Moreno, V. (2013) The mid p-value in exact tests for Hardy-Weinberg equilibrium, Statistical Applications in Genetics and Molecular Biology 12(4) pp. 433-448.

Graffelman J, Jain D, Weir B (2017) A genome-wide study of Hardy-Weinberg equilibrium with next generation sequence data. Human Genetics, 136, 727-741.

Examples

if (FALSE) { # \dontrun{
genome <- genometranslator::read_genome(
  data = "turtle.vcf",
  strata = "turtle.strata.tsv"
)

# Inspect HWE results and select thresholds interactively.
genome <- radr::filter_hwe(data = genome)

# Alternatively, start from a separate unfiltered GDS. Target markers in
# disequilibrium in four strata
# and use a mid-p threshold of 0.0001 (option 4).
scripted_genome <- genometranslator::read_genome("turtle_scripted.gds")
scripted_genome <- radr::filter_hwe(
  data = scripted_genome,
  interactive.filter = FALSE,
  filter.hwe = TRUE,
  hw.pop.threshold = 4,
  midp.threshold = 4L
)
} # }