Package {tidypopgen}

Title:	Tidy Population Genetics
Version:	1.0.0
Description:	We provide a tidy grammar of population genetics, facilitating the manipulation and analysis of data on biallelic single nucleotide polymorphisms (SNPs). 'tidypopgen' scales to very large genetic datasets by storing genotypes on disk, and performing operations on them in chunks, without ever loading all data in memory. The full functionalities of the package are described in Carter et al. (2025) <doi:10.1111/2041-210x.70204>.
License:	GPL (≥ 3)
Encoding:	UTF-8
Language:	en-GB
URL:	https://github.com/EvolEcolGroup/tidypopgen, https://evolecolgroup.github.io/tidypopgen/
BugReports:	https://github.com/EvolEcolGroup/tidypopgen/issues
Depends:	R (≥ 3.5.0), dplyr (≥ 1.2.0), tibble
Imports:	bigparallelr, bigsnpr, bigstatsr, data.table, foreach, generics, ggplot2 (≥ 3.5.0), methods, MASS, patchwork, runner, rlang, sf, stats, tidyselect, tidyr, utils, Rcpp, vctrs, gtable
Suggests:	adegenet, admixtools, broom, hexbin, hierfstat, knitr, detectRUNS, LEA, pcadapt, RhpcBLASctl, rmarkdown, rnaturalearth, rnaturalearthdata, readr, reticulate, testthat (≥ 3.0.0), vcfR, withr, xgboost, spelling
Additional_repositories:	https://evolecolgroup.r-universe.dev/
VignetteBuilder:	knitr
Config/testthat/edition:	3
LinkingTo:	Rcpp, RcppArmadillo (≥ 0.9.600), bigstatsr, rmio
LazyData:	true
Config/Needs/website:	rmarkdown
Config/roxygen2/version:	8.1.0
NeedsCompilation:	yes
Packaged:	2026-11-07 17:02:24 UTC; andrea
Author:	Evie Carter [aut], Eirlys Tysall [aut], Andrea Manica [aut, cre, cph], Chang Christopher [ctb] (Author of Hardy-Weinberg Equilibrium algorithm in PLINK 1.90, used in loci_hwe()), Shaun Purcell [ctb] (Author of Hardy-Weinberg Equilibrium algorithm in PLINK 1.90, used in loci_hwe()), Bengtsson Henrik [ctb] (Author of countLines in R.utils, modified for .vcf in count_vcf_variants())
Maintainer:	Andrea Manica <am315@cam.ac.uk>
Repository:	CRAN
Date/Publication:	2026-11-07 17:50:14 UTC

tidypopgen: Tidy Population Genetics

Description

We provide a tidy grammar of population genetics, facilitating the manipulation and analysis of data on biallelic single nucleotide polymorphisms (SNPs). 'tidypopgen' scales to very large genetic datasets by storing genotypes on disk, and performing operations on them in chunks, without ever loading all data in memory. The full functionalities of the package are described in Carter et al. (2025) doi:10.1111/2041-210x.70204.

Author(s)

Maintainer: Andrea Manica am315@cam.ac.uk (ORCID) [copyright holder]

Authors:

• Andrea Manica am315@cam.ac.uk (ORCID) [copyright holder]

• Evie Carter

• Eirlys Tysall

Other contributors:

• Chang Christopher (Author of Hardy-Weinberg Equilibrium algorithm in PLINK 1.90, used in loci_hwe()) [contributor]

• Shaun Purcell (Author of Hardy-Weinberg Equilibrium algorithm in PLINK 1.90, used in loci_hwe()) [contributor]

• Bengtsson Henrik (Author of countLines in R.utils, modified for .vcf in count_vcf_variants()) [contributor]

See Also

Useful links:

• https://github.com/EvolEcolGroup/tidypopgen

• https://evolecolgroup.github.io/tidypopgen/

• Report bugs at https://github.com/EvolEcolGroup/tidypopgen/issues

A $ method for gen_tibble objects

Description

A $ method for gen_tibble objects

Usage

## S3 replacement method for class 'gen_tbl'
x$i <- value

Arguments

Value

a gen_tibble

Examples

example_gt <- load_example_gt("gen_tbl")

# Add a new column
example_gt$region <- "East"

example_gt

Pipe operator

Description

See 'magrittr::pipe \

Usage

lhs %>% rhs

Arguments

Value

The result of calling rhs(lhs).

Examples

example_gt <- load_example_gt("gen_tbl")
example_gt %>% count_loci()

Subsetting method for matrix_int_names

Description

Subsetting method for matrix_int_names

Assignment method for matrix_int_names

Usage

## S3 method for class 'matrix_int_names'
x[i, j, i_names = NULL, j_names = NULL, drop = TRUE]

## S3 replacement method for class 'matrix_int_names'
x[i, j, i_names = NULL, j_names = NULL] <- value

Arguments

Value

A subsetted matrix_int_names object

See Also

Other matrix_int_names_functions: as.data.frame.matrix_int_names(), col_names(), matrix_int_names(), print.matrix_int_names(), row_names()

Examples

my_mat <- matrix_int_names(matrix(1:6, nrow = 2),
  row_names = c(10L, 20L),
  col_names = c(100L, 200L, 300L)
)
# subset by position
my_mat[1, 2:3, drop = FALSE]
# subset by integer names
my_mat[i_names = 20L, j_names = c(100L, 300L), drop = FALSE]
# assign by integer names
my_mat[i_names = 10L, j_names = 200L] <- 0L
my_mat

Attach a dendrogram panel to a heatmap

Description

Designed for use with |>. Accepts a ggplot or <heatmap_dendro> and returns a <heatmap_dendro> with the dendrogram panel inserted directly adjacent to the heatmap panel — pixel-perfect alignment regardless of axis label width. Can be chained to add dendrograms on both sides.

Usage

annotate_dendrogram(
  plot,
  hc,
  side = c("left", "top"),
  rel_size = 0.15,
  line_color = "grey30",
  line_size = 0.5
)

Arguments

Value

A <heatmap_dendro> object.

An arrange method for gen_tibble objects

Description

An arrange method for gen_tibble objects

Usage

## S3 method for class 'gen_tbl'
arrange(..., deparse.level = 1)

Arguments

Value

a gen_tibble

Examples

test_gt <- load_example_gt("gen_tbl")
test_gt %>% arrange(id)

An arrange method for grouped gen_tibble objects

Description

An arrange method for grouped gen_tibble objects

Usage

## S3 method for class 'grouped_gen_tbl'
arrange(..., deparse.level = 1)

Arguments

Value

a grouped gen_tibble

Examples

test_gt <- load_example_gt("grouped_gen_tbl")
test_gt %>% arrange(id)
test_gt <- load_example_gt("grouped_gen_tbl_sf")
test_gt %>% arrange(id)

Coercion method to data.frame for matrix_int_names

Description

This method converts a matrix_int_names object to a data.frame, preserving the integer row names. Otherwise, it works similarly to the default as.data.frame.matrix method.

Usage

## S3 method for class 'matrix_int_names'
as.data.frame(x, row.names = NULL, optional = FALSE, ...)

Arguments

Value

A data.frame with integer row names if present, otherwise character row names

See Also

Other matrix_int_names_functions: [.matrix_int_names(), col_names(), matrix_int_names(), print.matrix_int_names(), row_names()

Examples

# Create a matrix with integer row and column names
my_mat <- matrix_int_names(matrix(1:6, nrow = 2),
  row_names = c(10L, 20L),
  col_names = c(100L, 200L, 300L)
)
# Convert to data.frame
my_df <- as.data.frame(my_mat)
my_df

Convert runs of homozygosity to the format of detectRUNS

Description

Converts the runs found by windows_indiv_roh() to the data frame used by the detectRUNS package, so that its functions to plot and summarise runs (e.g. detectRUNS::plot_Runs(), detectRUNS::summaryRuns()) can be used.

Usage

as_detectRUNS(x)

Arguments

Value

a data frame with columns group (the grouping variables, pasted with "." if there are several, or the id if the runs were found on an ungrouped gen_tibble), id, chrom, nSNP, from, to and lengthBps, in the order of x

Examples

bed_file <- system.file("extdata", "roh", "roh_sim.bed",
  package = "tidypopgen"
)
sim_gt <- gen_tibble(bed_file, backingfile = tempfile(), quiet = TRUE)
sim_roh <- windows_indiv_roh(sim_gt %>% group_by(population))
head(as_detectRUNS(sim_roh))

Augment data with information from a gt_dapc object

Description

Augment for gt_dapc accepts a model object and a dataset and adds scores to each observation in the dataset. Scores for each component are stored in a separate column, which is given name with the pattern ".fittedLD1", ".fittedLD2", etc. For consistency with broom::augment.prcomp, a column ".rownames" is also returned; it is a copy of 'id', but it ensures that any scripts written for data augmented with broom::augment.prcomp will work out of the box (this is especially helpful when adapting plotting scripts).

Usage

## S3 method for class 'gt_dapc'
augment(x, data = NULL, k = NULL, ...)

Arguments

Value

A gen_tibble containing the original data along with additional columns containing each observation's projection into PCA space.

See Also

gt_dapc() gt_dapc_tidiers

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA and run DAPC
pca <- gt_pca_partialSVD(lobsters)
populations <- as.factor(lobsters$population)
dapc_res <- gt_dapc(pca, n_pca = 6, n_da = 2, pop = populations)

# Augment the gen_tibble with the DAPC scores
augment(dapc_res, data = lobsters, k = 2)

Augment data with information from a gt_pca object

Description

Augment for gt_pca accepts a model object and a dataset and adds scores to each observation in the dataset. Scores for each component are stored in a separate column, which is given name with the pattern ".fittedPC1", ".fittedPC2", etc. For consistency with broom::augment.prcomp, a column ".rownames" is also returned; it is a copy of 'id', but it ensures that any scripts written for data augmented with broom::augment.prcomp will work out of the box (this is especially helpful when adapting plotting scripts).

Usage

## S3 method for class 'gt_pca'
augment(x, data = NULL, k = NULL, ...)

Arguments

Value

A gen_tibble containing the original data along with additional columns containing each observation's projection into PCA space.

See Also

gt_pca_autoSVD() gt_pca_tidiers

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA object
pca <- gt_pca_partialSVD(lobsters)

# Augment the gen_tibble with PCA scores
augment(pca, data = lobsters)

# Adjust the number of components to add
augment(pca, data = lobsters, k = 2)

Augment the loci table with information from a analysis object

Description

augment_loci add columns to the loci table of a gen_tibble related to information from a given analysis.

Usage

augment_loci(x, data, ...)

Arguments

Value

A loci tibble with additional columns. If data is missing, a tibble of the information, with a column .rownames giving the loci names.

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA
pca <- gt_pca_partialSVD(lobsters)

# Augment the gen_tibble with the PCA scores
augment_loci(pca, data = lobsters)

Augment the loci table with information from a gt_pca object

Description

Augment for gt_pca accepts a model object and a gen_tibble and adds loadings for each locus to the loci table. Loadings for each component are stored in a separate column, which is given name with the pattern ".loadingPC1", ".loadingPC2", etc. If data is missing, then a tibble with the loadings is returned.

Usage

## S3 method for class 'gt_pca'
augment_loci(x, data = NULL, k = NULL, ...)

Arguments

Value

A gen_tibble with a loadings added to the loci tibble (accessible with show_loci(). If data is missing, a tibble of loadings.

See Also

gt_pca_autoSVD() gt_pca_tidiers

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA
pca <- gt_pca_partialSVD(lobsters)

# Augment the gen_tibble with the PCA scores
augment_loci(pca, data = lobsters)

Augment data with information from a q_matrix object

Description

Augment for q_matrix accepts a model object and a dataset and adds Q values to each observation in the dataset. Q values are stored in separate columns, which is given name with the pattern ".Q1",".Q2", etc. For consistency with broom::augment.prcomp, a column ".rownames" is also returned; it is a copy of 'id', but it ensures that any scripts written for data augmented with broom::augment.prcomp will work out of the box (this is especially helpful when adapting plotting scripts).

Usage

## S3 method for class 'q_matrix'
augment(x, data = NULL, ...)

Arguments

Value

A gen_tibble containing the original data along with additional columns containing each observation's Q values.

Examples

# run the example only if we have the package installed
if (requireNamespace("LEA", quietly = TRUE)) {
  example_gt <- load_example_gt("gen_tbl")

  # Create a gt_admix object
  admix_obj <- example_gt %>% gt_snmf(k = 1:3, project = "force")

  # Extract a Q matrix
  q_mat_k3 <- get_q_matrix(admix_obj, k = 3, run = 1)

  # Augment the gen_tibble with Q values
  augment(q_mat_k3, data = example_gt)
}

Autoplots for gt_cluster_pca objects

Description

For gt_cluster_pca, autoplot produces a plot of a metric of choice ('BIC', 'AIC' or 'WSS') against the number of clusters (k). This plot is can be used to infer the best value of k, which corresponds to the smallest value of the metric (the minimum in an 'elbow' shaped curve). In some cases, there is not 'elbow' and the metric keeps decreasing with increasing k; in such cases, it is customary to choose the value of k at which the decrease in the metric reaches as plateau. For a programmatic way of choosing k, use gt_cluster_pca_best_k().

Usage

## S3 method for class 'gt_cluster_pca'
autoplot(object, metric = c("BIC", "AIC", "WSS"), ...)

Arguments

Details

autoplot produces simple plots to quickly inspect an object. They are not customisable; we recommend that you use ggplot2 to produce publication ready plots.

Value

a ggplot2 object

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA object
pca <- gt_pca_partialSVD(lobsters)

# Run clustering on the first 10 PCs
cluster_pca <- gt_cluster_pca(
  x = pca,
  n_pca = 10,
  k_clusters = c(1, 5),
  method = "kmeans",
  n_iter = 1e5,
  n_start = 10,
  quiet = FALSE
)

# Autoplot BIC
autoplot(cluster_pca, metric = "BIC")

# # Autoplot AIC
autoplot(cluster_pca, metric = "AIC")

# # Autoplot WSS
autoplot(cluster_pca, metric = "WSS")

Autoplots for gt_dapc objects

Description

For gt_dapc, the following types of plots are available:

• screeplot: a plot of the eigenvalues of the discriminant axes

• scores a scatterplot of the scores of each individual on two discriminant axes (defined by ld)

• loadings a plot of loadings of all loci for a discriminant axis (chosen with ld)

• components a bar plot showing the probability of assignment to each cluster

Usage

## S3 method for class 'gt_dapc'
autoplot(
  object,
  type = c("screeplot", "scores", "loadings", "components"),
  ld = NULL,
  group = NULL,
  n_col = 1,
  ...
)

Arguments

Details

autoplot produces simple plots to quickly inspect an object. They are not customisable; we recommend that you use ggplot2 to produce publication ready plots.

Value

a ggplot2 object

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA and run DAPC
pca <- gt_pca_partialSVD(lobsters)
populations <- as.factor(lobsters$population)
dapc_res <- gt_dapc(pca, n_pca = 6, n_da = 2, pop = populations)

# Screeplot
autoplot(dapc_res, type = "screeplot")

# Scores plot
autoplot(dapc_res, type = "scores", ld = c(1, 2))

# Loadings plot
autoplot(dapc_res, type = "loadings", ld = 1)

# Components plot
autoplot(dapc_res, type = "components", group = populations)

Autoplots for indiv_roh objects

Description

Plots of runs of homozygosity found by windows_indiv_roh(), equivalent to the plots of detectRUNS. Runs are coloured (or split) by group if the gen_tibble was grouped. The following types of plots are available:

• runs: the runs of each individual along each chromosome (as detectRUNS::plot_Runs());

• stacked: the runs of each group stacked along each chromosome, one run per line, ordered by start (as detectRUNS::plot_StackedRuns());

• coverage: the percentage of individuals of each group in a run along each chromosome (as detectRUNS::plot_SnpsInRuns() and detectRUNS::plot_manhattanRuns(), computed from the runs rather than at each locus);

• froh: the inbreeding coefficient F_ROH, the total length of the runs of an individual over the length of the genome (with froh_by = "individual", one point per individual, by group) or of each chromosome (with froh_by = "chromosome", the mean of each group) (as detectRUNS::plot_InbreedingChr()); the length of a chromosome is its largest position, as in detectRUNS;

• violin: the distribution of the total (or mean, see method) length of the runs of each individual, by group (as detectRUNS::plot_ViolinRuns());

• length_class: the percentage of runs in each class of length, by group (as detectRUNS::plot_DistributionRuns()).

Usage

## S3 method for class 'indiv_roh'
autoplot(
  object,
  type = c("runs", "stacked", "coverage", "froh", "violin", "length_class"),
  chromosome = NULL,
  froh_by = c("individual", "chromosome"),
  method = c("sum", "mean"),
  classes = c(0, 2, 4, 8, 16),
  ...
)

Arguments

Details

Individuals with no runs count in coverage, froh and violin (with no runs, or a length of 0).

autoplot produces simple plots to quickly inspect an object. They are not customisable; we recommend that you use ggplot2 to produce publication ready plots.

Value

a ggplot2 object

Examples

# a simulated dataset of three populations with increasing inbreeding
bed_file <- system.file("extdata", "roh", "roh_sim.bed",
  package = "tidypopgen"
)
sim_gt <- gen_tibble(bed_file, backingfile = tempfile(), quiet = TRUE)
sim_roh <- windows_indiv_roh(sim_gt %>% group_by(population))

autoplot(sim_roh, type = "runs", chromosome = "3")
autoplot(sim_roh, type = "coverage")
autoplot(sim_roh, type = "froh")
autoplot(sim_roh, type = "length_class")

Autoplots for qc_report_indiv objects

Description

For qc_report_indiv, the following types of plots are available:

• scatter: a plot of missingness and observed heterozygosity within individuals.

• relatedness: a histogram of paired kinship coefficients

• histogram: for gen_tibbles containing pseudohaploid data, a histogram of missingness, split by ploidy.

Usage

## S3 method for class 'qc_report_indiv'
autoplot(
  object,
  type = c("scatter", "relatedness", "histogram"),
  miss_threshold = 0.05,
  kings_threshold = NULL,
  ...
)

Arguments

Details

autoplot produces simple plots to quickly inspect an object. They are not customisable; we recommend that you use ggplot2 to produce publication ready plots.

Value

a ggplot2 object. For relatedness plots of a report created from a grouped gen_tibble, a named list of ggplot2 objects, one for each group.

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
example_gt <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Create QC report for individuals
indiv_report <- example_gt %>% qc_report_indiv()

# Autoplot missingness and observed heterozygosity
autoplot(indiv_report, type = "scatter", miss_threshold = 0.1)

# Create QC report with kinship filtering
indiv_report_rel <-
  example_gt %>% qc_report_indiv(kings_threshold = "second")

# Autoplot relatedness
autoplot(indiv_report_rel, type = "relatedness", kings_threshold = "second")

Autoplots for qc_report_loci objects

Description

For qc_report_loci, the following types of plots are available:

• overview: an UpSet plot, giving counts of snps over the threshold for missingness, minor allele frequency, and Hardy-Weinberg equilibrium P-value, and visualising the interaction between these.

• all: a four panel plot, containing ⁠missing high maf⁠, ⁠missing low maf⁠, hwe, and ⁠significant hwe⁠ plots

• missing: a histogram of proportion of missing data

• ⁠missing low maf⁠: a histogram of the proportion of missing data for snps with low minor allele frequency

• ⁠missing high maf⁠:a histogram of the proportion of missing data for snps with high minor allele frequency

• maf: a histogram of minor allele frequency

• hwe: a histogram of HWE exact test p-values

• ⁠significant hwe⁠: a histogram of significant HWE exact test p-values

Usage

## S3 method for class 'qc_report_loci'
autoplot(
  object,
  type = c("overview", "all", "missing", "missing low maf", "missing high maf", "maf",
    "hwe", "significant hwe"),
  maf_threshold = 0.05,
  miss_threshold = 0.01,
  hwe_p = 0.01,
  ...
)

Arguments

Details

autoplot produces simple plots to quickly inspect an object. They are not customisable; we recommend that you use ggplot2 to produce publication ready plots.

Value

a ggplot2 object

Examples

# Create a gen_tibble
bed_file <-
  system.file("extdata", "related", "families.bed", package = "tidypopgen")
example_gt <- gen_tibble(bed_file,
  backingfile = tempfile("families"),
  quiet = TRUE,
  valid_alleles = c("1", "2")
)

loci_report <- example_gt %>% qc_report_loci()

# Plot the QC report overview
autoplot(loci_report, type = "overview")

# Plot the QC report all
autoplot(loci_report, type = "all")

# Plot missing data
autoplot(loci_report, type = "missing")

# Plot missing with low maf
autoplot(loci_report, type = "missing low maf", maf_threshold = 0.05)

# Plot missing with high maf
autoplot(loci_report, type = "missing high maf", maf_threshold = 0.05)

# Plot maf
autoplot(loci_report, type = "maf", maf_threshold = 0.05)

# Plot hwe
autoplot(loci_report, type = "hwe", hwe_p = 0.01)

# Plot significant hwe
autoplot(loci_report, type = "significant hwe", hwe_p = 0.01)

Autoplots for gt_admix objects

Description

For gt_admix, the following types of plots are available:

• cv: the cross-validation error for each value of k

• barplot a standard barplot of the admixture proportions

Usage

## S3 method for class 'gt_admix'
autoplot(object, type = c("cv", "barplot"), k = NULL, run = NULL, ...)

Arguments

Details

autoplot produces simple plots to quickly inspect an object. They are not customisable; we recommend that you use ggplot2 to produce publication ready plots.

This autoplot will automatically rearrange individuals according to their id and any grouping variables if an associated 'data' gen_tibble is provided. To avoid any automatic re-sorting of individuals, set arrange_by_group and arrange_by_indiv to FALSE. See autoplot.q_matrix for further details.

Value

a ggplot2 object

Examples

# Read example gt_admix object
admix_obj <-
  readRDS(system.file("extdata", "anolis", "anole_adm_k3.rds",
    package = "tidypopgen"
  ))
# Cross-validation plot
autoplot(admix_obj, type = "cv")

# Basic barplot
autoplot(admix_obj, k = 3, run = 1, type = "barplot")

# Barplot with individuals arranged by Q proportion
# (using additional arguments, see `autoplot.q_matrix` for details)
autoplot(admix_obj,
  k = 3, run = 1, type = "barplot", annotate_group = TRUE,
  arrange_by_group = TRUE, arrange_by_indiv = TRUE,
  reorder_within_groups = TRUE
)

Autoplots for gt_pca objects

Description

For gt_pca, the following types of plots are available:

• screeplot: a plot of the eigenvalues of the principal components (currently it plots the singular value)

• scores a scatterplot of the scores of each individual on two principal components (defined by k)

• loadings a plot of loadings of all loci for a given component (chosen with k)

Usage

## S3 method for class 'gt_pca'
autoplot(object, type = c("screeplot", "scores", "loadings"), k = NULL, ...)

Arguments

Details

autoplot produces simple plots to quickly inspect an object. They are not customisable; we recommend that you use ggplot2 to produce publication ready plots.

Value

a ggplot2 object

Examples

library(ggplot2)
# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA object
pca <- gt_pca_partialSVD(lobsters)

# Screeplot
autoplot(pca, type = "screeplot")

# Scores plot
autoplot(pca, type = "scores")

# Colour by population
autoplot(pca, type = "scores") + aes(colour = lobsters$population)

# Scores plot of different components
autoplot(pca, type = "scores", k = c(1, 3)) +
  aes(colour = lobsters$population)

Autoplots for gt_pcadapt objects

Description

For gt_pcadapt, the following types of plots are available:

• qq: a quantile-quantile plot of the p-values from pcadapt (wrapping bigsnpr::snp_qq())

• manhattan a manhattan plot of the p-values from pcadapt (wrapping bigsnpr::snp_manhattan())

Usage

## S3 method for class 'gt_pcadapt'
autoplot(object, type = c("qq", "manhattan"), ...)

Arguments

Details

autoplot produces simple plots to quickly inspect an object. They are not customisable; we recommend that you use ggplot2 to produce publication ready plots.

Value

a ggplot2 object

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA object
pca <- gt_pca_partialSVD(lobsters)

# Create a gt_pcadapt object
pcadapt_obj <- gt_pcadapt(lobsters, pca, k = 2)

# Plot the p-values from pcadapt
autoplot(pcadapt_obj, type = "qq")

# Plot the manhattan plot of the p-values from pcadapt
autoplot(pcadapt_obj, type = "manhattan")

Autoplot for pairwise distance matrices

Description

This function provides a convenient way to visualize pairwise distance matrices (e.g., Fst, IBS) as heatmaps using ggplot2. It supports both matrix and tidy tibble formats, and can optionally reorder rows/columns using either an explicit ordering vector or a function (e.g. hierarchical clustering).

Usage

## S3 method for class 'pairwise_matrix'
autoplot(object, order = NULL, ...)

## S3 method for class 'pairwise_tbl'
autoplot(object, order = NULL, ...)

Arguments

Value

A ggplot object.

Examples

example_gt <- load_example_gt("gen_tbl")

# With a pairwise matrix:
fst_dist <- example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Nei87", type = "pairwise")

autoplot(fst_dist)

# Explicit ordering
autoplot(
  fst_dist,
  order = rev(rownames(fst_dist))
)

# Hierarchical clustering
autoplot(
  fst_dist,
  order = function(x) {
    stats::hclust(stats::as.dist(x), method = "complete")
  }
)

# With a tidy tibble of distances:
fst_tidy <- example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Nei87", type = "tidy")

autoplot(fst_tidy)

autoplot(
  fst_tidy,
  order = function(x) {
    stats::hclust(stats::as.dist(x))
  }
)

Autoplots for q_matrix objects

Description

This autoplot will automatically rearrange individuals according to their id and any grouping variables if an associated 'data' gen_tibble is provided. To avoid any automatic re-sorting of individuals, set arrange_by_group and arrange_by_indiv to FALSE.

Usage

## S3 method for class 'q_matrix'
autoplot(
  object,
  data = NULL,
  annotate_group = TRUE,
  arrange_by_group = TRUE,
  arrange_by_indiv = TRUE,
  reorder_within_groups = FALSE,
  ...
)

Arguments

Value

a barplot of individuals, coloured by ancestry proportion

Examples

# Read example gt_admix obejct
admix_obj <-
  readRDS(system.file("extdata", "anolis", "anole_adm_k3.rds",
    package = "tidypopgen"
  ))

# Extract a Q matrix
q_mat_k3 <- get_q_matrix(admix_obj, k = 3, run = 1)

# Basic autoplot
autoplot(q_mat_k3, annotate_group = FALSE, arrange_by_group = FALSE)

# To arrange individuals by group and by Q proportion
autoplot(q_mat_k3,
  annotate_group = TRUE, arrange_by_group = TRUE,
  arrange_by_indiv = TRUE, reorder_within_groups = TRUE
)

Plot site frequency spectra

Description

Plot the site frequency spectra computed by pop_sfs() (bar plots of the number of sites by allele count, one panel per population), pairwise_pop_sfs() (heatmaps of the joint spectra, one panel per pair of populations), or a single spectrum returned with type = "array" (a bar plot for one population, a heatmap for two). The output of pop_sfs_preview() is plotted as the number of segregating sites for each projection, one panel per population, with the projection chosen by default marked. Spectra of more than two populations (from nwise_pop_sfs()) cannot be plotted; use pairwise_pop_sfs() to plot the joint spectra of pairs of populations.

Usage

## S3 method for class 'sfs'
autoplot(object, masked = FALSE, log = NULL, min_sites = NULL, ...)

## S3 method for class 'pop_sfs_tbl'
autoplot(object, masked = FALSE, log = NULL, ...)

## S3 method for class 'pairwise_sfs_tbl'
autoplot(object, masked = FALSE, log = NULL, min_sites = NULL, ...)

## S3 method for class 'sfs_preview_tbl'
autoplot(object, ...)

Arguments

Details

The masked cells of the spectra (the corners, i.e. sites monomorphic in the projected sample, and the cells that were folded out in a folded spectrum; see pop_sfs()) are not shown by default. In heatmaps, the numbers of sites are shown on a log scale by default, and empty cells are left blank. As the projection spreads each locus over many cells, many cells have a tiny fraction of a site: on a log scale, the colour scale starts at min_sites (by default, a thousandth of the largest number of sites in a cell), and cells with fewer sites are shown with its lightest colour (as vmin in the plots of dadi).

Value

a ggplot2 object.

See Also

pop_sfs(), pairwise_pop_sfs(), pop_sfs_preview()

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# the 1D spectrum of each population
example_gt %>%
  pop_sfs(n_proj = 4) %>%
  autoplot()

# the joint spectrum of each pair of populations
example_gt %>%
  pairwise_pop_sfs(n_proj = 4, folded = FALSE) %>%
  autoplot()

# a single spectrum
sfs_2d <- example_gt %>%
  pairwise_pop_sfs(n_proj = 4, folded = FALSE, type = "array")
autoplot(sfs_2d$pop1.pop2)

# the number of segregating sites for each projection
example_gt %>%
  pop_sfs_preview() %>%
  autoplot()

Block bootstrap of site frequency spectra

Description

Compute bootstrap replicates of the site frequency spectra of pop_sfs(), pairwise_pop_sfs() and nwise_pop_sfs(), by resampling blocks of loci with replacement. Replicates are used to estimate the uncertainty of statistics computed from the spectra, or of the parameters of demographic models fitted with dadi, moments or fastsimcoal2 (see write_sfs() to write the replicates in their formats).

Usage

bootstrap_pop_sfs(.x, ...)

## S3 method for class 'gen_tbl'
bootstrap_pop_sfs(
  .x,
  n_bootstrap = 100,
  block_size,
  size_unit = c("snp", "cM", "bp"),
  n_proj = NULL,
  folded = TRUE,
  ancestral = NULL,
  type = c("tidy", "array"),
  n_cores = bigstatsr::nb_cores(),
  ...
)

bootstrap_pairwise_pop_sfs(.x, ...)

## S3 method for class 'gen_tbl'
bootstrap_pairwise_pop_sfs(
  .x,
  n_bootstrap = 100,
  block_size,
  size_unit = c("snp", "cM", "bp"),
  n_proj = NULL,
  folded = TRUE,
  ancestral = NULL,
  type = c("tidy", "array"),
  fastest = c("last", "first"),
  n_cores = bigstatsr::nb_cores(),
  ...
)

bootstrap_nwise_pop_sfs(.x, ...)

## S3 method for class 'gen_tbl'
bootstrap_nwise_pop_sfs(
  .x,
  n_bootstrap = 100,
  block_size,
  size_unit = c("snp", "cM", "bp"),
  n_proj = NULL,
  folded = TRUE,
  ancestral = NULL,
  type = c("tidy", "array"),
  fastest = c("last", "first"),
  max_cells = 1e+07,
  n_cores = bigstatsr::nb_cores(),
  ...
)

Arguments

Details

The loci are divided into non-overlapping blocks that never span more than one chromosome, defined as in the ⁠jackknife_*⁠ functions (e.g. jackknife_pop_tajimas_d()): blocks of block_size SNPs, base pairs or centiMorgans (the latter requires the genetic_dist column of the loci table). Blocks should be larger than the extent of linkage disequilibrium, so that they are approximately independent. Each replicate draws, with replacement, as many blocks as there are (as the bootstrap of dadi, Misc.bootstraps_from_dd_chunks()). With block_size = 1 and size_unit = "snp", single SNPs are resampled (as the --resample-vcf option of easySFS), which ignores linkage.

A spectrum is a sum over loci, so the spectrum of a replicate is the sum of the spectra of the blocks it drew (counting each block as many times as it was drawn). The genotypes are therefore read only once, to compute the spectrum of each block. Missing data are dealt with by projection, as in pop_sfs(), with the same projections (n_proj) for all replicates: if n_proj is NULL, they are chosen from all the loci with pop_sfs_preview(). The blocks are drawn with the random number generator of R, so use set.seed() to make the replicates reproducible.

Value

If type = "tidy", a tibble with the columns of the tidy output of pop_sfs(), pairwise_pop_sfs() or nwise_pop_sfs(), and a column replicate (after the group columns) with the number of the replicate. If type = "array", a list with one element per replicate, each of which is what pop_sfs(), pairwise_pop_sfs() or nwise_pop_sfs() return with type = "array" (lists of sfs objects, or an sfs object). In both cases, the attribute bootstrap_blocks is a matrix with the number of times each block (columns, numbered in the order in which they appear in the loci table) was drawn in each replicate (rows), and the tibble also has the attributes n_proj and folded.

See Also

pop_sfs(), write_sfs(), jackknife_stats_generic()

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

set.seed(1)
# resample blocks of 2 SNPs
boot <- example_gt %>%
  bootstrap_pop_sfs(n_bootstrap = 10, block_size = 2, n_proj = 2)
boot

# the 95% interval of the number of sites in each cell
boot %>%
  dplyr::group_by(group, allele_count) %>%
  dplyr::summarise(
    lower = stats::quantile(n_sites, 0.025),
    upper = stats::quantile(n_sites, 0.975)
  )

# the joint spectra of each pair, as arrays
boot_2d <- example_gt %>%
  bootstrap_pairwise_pop_sfs(
    n_bootstrap = 5, block_size = 2, n_proj = 2,
    type = "array"
  )
boot_2d[[1]]$pop1.pop2

Combine method for gt_admix objects

Description

Combine method for gt_admix objects

Usage

## S3 method for class 'gt_admix'
c(..., match_attributes = TRUE)

Arguments

Value

A gt_admix object with the combined data

Examples

# run the example only if we have the package installed
if (requireNamespace("LEA", quietly = TRUE)) {
  example_gt <- load_example_gt("gen_tbl")

  # Create a gt_admix object
  admix_obj <- example_gt %>% gt_snmf(k = 1:3, project = "force")

  # Create a second gt_admix object
  admix_obj2 <- example_gt %>% gt_snmf(k = 2:4, project = "force")

  # Combine the two gt_admix objects
  new_admix_obj <- c(admix_obj, admix_obj2)
  summary(new_admix_obj)
}

Combine a gen_tibble to a data.frame or tibble by column

Description

A cbind() method to merge gen_tibble objects with data.frames and normal tibbles. Whilst this works, it is not ideal as it does not check the order of the tables, and we suggest that you use dplyr::left_join() instead. Note that cbind will not combine two gen_tibbles (i.e. it will NOT combine markers for the same individuals)

Usage

## S3 method for class 'gen_tbl'
cbind(..., deparse.level = 1)

## S3 method for class 'grouped_gen_tbl'
cbind(..., deparse.level = 1)

Arguments

Details

As in dplyr::bind_cols() (which dplyr uses to cbind() grouped tibbles), duplicated column names are made unique (e.g. population...2 and population...4), with a message. If a grouping variable is renamed, the result is no longer grouped. The columns id and genotypes of the gen_tibble can not be duplicated, as renaming them would break it.

Value

a gen_tibble

Examples

example_gt <- load_example_gt("gen_tbl")

# Create a dataframe to combine with the gen_tibble
df <- data.frame(region = c("A", "A", "B", "B", "A", "B", "B"))

# Combine the gen_tibble with the dataframe
example_gt <- cbind(example_gt, df)

Get or set column names for matrix_int_names

Description

This function reliably returns integer names if present, otherwise character dimnames. Use this instead of colnames() if you need guaranteed dispatch.

Usage

col_names(x)

## Default S3 method:
col_names(x)

## S3 method for class 'matrix_int_names'
col_names(x)

col_names(x) <- value

## Default S3 replacement method:
col_names(x) <- value

## S3 replacement method for class 'matrix_int_names'
col_names(x) <- value

Arguments

Value

Integer or character vector of column names, or NULL

See Also

Other matrix_int_names_functions: [.matrix_int_names(), as.data.frame.matrix_int_names(), matrix_int_names(), print.matrix_int_names(), row_names()

Examples

my_mat <- matrix_int_names(matrix(1:6, nrow = 2),
  row_names = c(10L, 20L),
  col_names = c(100L, 200L, 300L)
)
col_names(my_mat)
# set character names instead
col_names(my_mat) <- c("a", "b", "c")
col_names(my_mat)

Count the number of loci in a gen_tibble

Description

Count the number of loci in gen_tibble (or directly from its genotype column).

Usage

count_loci(.x, ...)

## S3 method for class 'tbl_df'
count_loci(.x, ...)

## S3 method for class 'vctrs_bigSNP'
count_loci(.x, ...)

Arguments

Value

the number of loci

Examples

example_gt <- load_example_gt("gen_tbl")

example_gt %>% count_loci()

Distruct colours

Description

Colours in the palette used by distruct

Usage

distruct_colours

Format

A vector of 56 hex colours

Details

The 56 distinct colours of the palette of distruct (which repeats 4 of its 60 colours), in their original order. scale_fill_distruct() recycles them when there are more groups.

Tidyverse methods for gt objects

Description

A filter method for gen_tibble objects

Usage

## S3 method for class 'gen_tbl'
filter(..., deparse.level = 1)

Arguments

Value

a gen_tibble

Examples

test_gt <- load_example_gt("gen_tbl")
test_gt %>% filter(id %in% c("a", "c"))

A filter method for grouped gen_tibble objects

Description

A filter method for grouped gen_tibble objects

Usage

## S3 method for class 'grouped_gen_tbl'
filter(..., deparse.level = 1)

Arguments

Value

a grouped gen_tibble

Examples

test_gt <- load_example_gt("grouped_gen_tbl")
test_gt %>% filter(id %in% c("a", "c"))
test_gt <- load_example_gt("grouped_gen_tbl_sf")
test_gt %>% filter(id %in% c("a", "c"))

Filter individuals based on a relationship threshold

Description

This function takes a matrix of x by y individuals containing relatedness coefficients and returns the maximum set of individuals that contains no relationships above the given threshold.

Usage

filter_high_relatedness(
  matrix,
  .x = NULL,
  kings_threshold = NULL,
  verbose = FALSE
)

Arguments

Details

Individuals are removed greedily, as in caret::findCorrelation(): individuals are considered in order of decreasing mean relatedness, and for each pair related above kings_threshold, the first individual is removed if its mean relatedness is larger than the mean of the remaining matrix without the second individual, and the second one otherwise. No pair above the threshold is left, but the set of individuals kept is not guaranteed to be the largest possible.

Missing relationship coefficients (e.g. KING kinship for pairs of individuals with no loci genotyped in both) are not evidence of relatedness: such pairs are treated as not related, with a warning.

Value

a list where '1' is individual ID's to retain, '2' is individual ID's to remove, and '3' is a boolean where individuals to keep are TRUE and individuals to remove are FALSE

Examples

example_gt <- load_example_gt("gen_tbl")

# Calculate relationship matrix
king_matrix <- example_gt %>% pairwise_king(as_matrix = TRUE)

# Filter individuals with threshold above 0.2
filter_high_relatedness(king_matrix, example_gt, kings_threshold = 0.2)

Find duplicates in the loci table

Description

This function finds duplicated SNPs by checking the positions within each chromosome. It can return a list of duplicated SNPs or a logical value indicating whether there are any duplicated loci.

Usage

find_duplicated_loci(.x, error_on_false = FALSE, list_duplicates = TRUE, ...)

Arguments

Value

If list_duplicates is TRUE, returns a character vector of duplicated loci names (character(0) when none). If list_duplicates is FALSE, returns TRUE when no duplicates exist and FALSE when duplicates are present. If error_on_false is TRUE and duplicates exist, an error is thrown.

Examples

example_gt <- load_example_gt("gen_tbl")
show_loci(example_gt) <- test_loci <- data.frame(
  big_index = c(1:6),
  name = paste0("rs", 1:6),
  chromosome = paste0("chr", c(1, 1, 1, 1, 1, 1)),
  position = as.integer(c(3, 3, 5, 65, 343, 46)),
  genetic_dist = as.double(rep(0, 6)),
  allele_ref = c("A", "T", "C", "G", "C", "T"),
  allele_alt = c("T", "C", NA, "C", "G", "A")
)

show_loci(example_gt)

# Find which loci are duplicated
example_gt %>% find_duplicated_loci()

Constructor for a gen_tibble

Description

A gen_tibble stores genotypes for individuals in a tidy format. DESCRIBE here the format

Usage

gen_tibble(
  x,
  ...,
  valid_alleles = c("A", "T", "C", "G"),
  missing_alleles = c("0", "."),
  backingfile = NULL,
  allow_duplicates = FALSE,
  quiet = FALSE
)

## S3 method for class 'character'
gen_tibble(
  x,
  ...,
  parser = c("cpp", "vcfR"),
  n_cores = 1,
  chunk_size = NULL,
  valid_alleles = c("A", "T", "C", "G"),
  missing_alleles = c("0", "."),
  backingfile = NULL,
  allow_duplicates = FALSE,
  names_as_int = FALSE,
  quiet = FALSE
)

## S3 method for class 'matrix'
gen_tibble(
  x,
  indiv_meta,
  loci,
  ...,
  ploidy = 2,
  valid_alleles = c("A", "T", "C", "G"),
  missing_alleles = c("0", "."),
  backingfile = NULL,
  allow_duplicates = FALSE,
  quiet = FALSE
)

Arguments

Details

• VCF files: the fast cpp parser is used by default. Both cpp and vcfR parsers attempt to establish ploidy from the first variant; if that variant is found in a sex chromosome (or mtDNA), the parser will fail with 'Error: a genotype has more than max_ploidy alleles...'. To successful import such a VCF, change the order of variants so that the first chromosome is an autosome using a tool such as vcftools. Individuals with a missing genotype ('.') on the first variant get their ploidy from their first genotype that is not '.': the cpp parser looks for it in the next 1,000 variants only (so that the search does not slow down parsing), the vcfR parser in the whole file. Individuals with no such genotype get the largest ploidy of the other individuals (or 2, if no ploidy is known). Currently, only biallelic SNPs are supported. If haploid variants (e.g. sex chromosomes) are included in the VCF, they are not transformed into homozygous calls. Instead, reference alleles will be coded as 0 and alternative alleles will be coded as 1.

• PLINK files (.bed, .ped, or a bigSNP .rds): the individual information (from the .fam file, or the first six columns of the .ped file) is stored in the same columns for all formats: the individual ID as id; the family ID as population (unless it is the same as the individual ID for all individuals); the paternal and maternal IDs as paternal_ID and maternal_ID; sex as sex, a factor with levels "female" (PLINK code 2) and "male" (code 1); and the phenotype as phenotype, a factor with levels "control" (code 1) and "case" (code 2), or a number for a quantitative phenotype (with values other than 1, 2, 0 and -9). Missing values (0 for parents and sex, 0 or -9 for a case/control phenotype, and -9 for a quantitative one) are NA, and paternal_ID, maternal_ID, sex and phenotype are only created if they have some non-missing values. gt_as_plink() writes these columns back to PLINK files.

• packedancestry files: When loading packedancestry files, missing alleles will be converted from 'X' to NA

Value

an object of the class gen_tbl.

Note

Helper functions for accessing gen_tibble object attributes and checking gen_tibble ploidy can be found in gt_helper_functions.R

Examples

# Create a gen_tibble from a .bed file
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Create a gen_tibble from a .vcf file
vcf_path <-
  system.file("extdata", "anolis",
    "punctatus_t70_s10_n46_filtered.recode.vcf.gz",
    package = "tidypopgen"
  )
gen_tibble(vcf_path, quiet = TRUE, backingfile = tempfile("anolis_"))

# Create a gen_tibble from a matrix of genotypes:
test_indiv_meta <- data.frame(
  id = c("a", "b", "c"),
  population = c("pop1", "pop1", "pop2")
)
test_genotypes <- rbind(
  c(1, 1, 0, 1, 1, 0),
  c(2, 1, 0, 0, 0, 0),
  c(2, 2, 0, 0, 1, 1)
)
test_loci <- data.frame(
  name = paste0("rs", 1:6),
  chromosome = paste0("chr", c(1, 1, 1, 1, 2, 2)),
  position = as.integer(c(3, 5, 65, 343, 23, 456)),
  genetic_dist = as.double(rep(0, 6)),
  allele_ref = c("A", "T", "C", "G", "C", "T"),
  allele_alt = c("T", "C", NA, "C", "G", "A")
)

gen_tibble(
  x = test_genotypes,
  loci = test_loci,
  indiv_meta = test_indiv_meta,
  valid_alleles = c("A", "T", "C", "G"),
  quiet = TRUE
)

Restore the class of a tibble from its genotypes column

Description

Some dplyr operations return tibbles that have lost the class of a gen_tibble, such as the data frames returned by pick() or across() within summarise(), reframe() or mutate(), the result of reframe(), or as_tibble(). Their genotypes column, however, keeps its class. The generics of tidypopgen therefore have fallback methods for data frames (and grouped data frames), which use restore_genotypes_tbl() to give back the tibble its class, as given by genotypes_tbl_class() for its genotypes column, and then dispatch again.

Usage

genotypes_tbl_class(genotypes, ...)

## Default S3 method:
genotypes_tbl_class(genotypes, ...)

## S3 method for class 'vctrs_bigSNP'
genotypes_tbl_class(genotypes, ...)

restore_genotypes_tbl(.x, generic = "this function")

Arguments

Details

Packages that build their own tibbles on tidypopgen (e.g. the gl_tibble of tidypopgl) add a method to genotypes_tbl_class() for the class of their genotypes column, so that their own methods of the tidypopgen generics are used after the class is restored.

Value

genotypes_tbl_class() returns the class of the tibble that holds the genotypes column ("gen_tbl" for a vctrs_bigSNP column), or NULL for a column of an unknown class.

restore_genotypes_tbl() returns .x as a tibble with the class given by genotypes_tbl_class(), and its grouped version (e.g. "grouped_gen_tbl") if .x is grouped. A rowwise_df is ungrouped. It stops with an error if .x already has that class, as the generic would otherwise dispatch to the fallback method again.

Return a single P matrix from a gt_admix object

Description

This function retrieves a single P matrix from a gt_admix object based on the specified k value and run number.

Usage

get_p_matrix(x, ..., k, run)

Arguments

Value

A single P matrix from the gt_admix object

Examples

# Read example gt_admix object
admix_obj <-
  readRDS(system.file("extdata", "anolis", "anole_adm_k3.rds",
    package = "tidypopgen"
  ))

# Extract a P matrix
get_p_matrix(admix_obj, k = 3, run = 1)

Return a single Q matrix from a gt_admix object

Description

This function retrieves a single Q matrix from a gt_admix object based on the specified k value and run number.

Usage

get_q_matrix(x, ..., k, run)

Arguments

Value

A single Q matrix from the gt_admix object

Examples

# Read example gt_admix obejct
admix_obj <-
  readRDS(system.file("extdata", "anolis", "anole_adm_k3.rds",
    package = "tidypopgen"
  ))

# Extract a Q matrix
get_q_matrix(admix_obj, k = 3, run = 1)

Add an simple feature geometry to a gen_tibble

Description

gt_add_sf adds an active sf geometry column to a gen_tibble object. The resulting gen_tbl inherits from sf and can be used with functions from the sf package. It is possible to either create a sf::sfc geometry column from coordinates, or to provide an existing geometry column (which will then become the active geometry for sf).

Usage

gt_add_sf(x, coords = NULL, crs = NULL, sfc_column = NULL)

Arguments

Value

a gen_tibble object with an additional geometry column (and thus belonging also to sf class).

Examples

example_gt <- load_example_gt("gen_tbl")

# Add some coordinates
example_gt <- example_gt %>% mutate(
  longitude = c(0, 0, 2, 2, 0, 2, 2),
  latitude = c(51, 51, 49, 49, 51, 41, 41)
)

# Convert lat and long to sf:
example_gt <- gt_add_sf(x = example_gt, coords = c("longitude", "latitude"))

# Check class
class(example_gt)

Reorder the q matrices based on the grouping variable

Description

This function reorders the q matrices in a gt_admix object based on the grouping variable. This is useful before plotting when the samples from each group are not adjacent to each other in the q matrix.

Usage

gt_admix_reorder_q(x, group = NULL)

Arguments

Value

a gt_admix object with the q matrices reordered

Examples

# run the example only if we have the package installed
if (requireNamespace("LEA", quietly = TRUE)) {
  example_gt <- load_example_gt("gen_tbl")

  # Create a gt_admix object
  admix_obj <- example_gt %>% gt_snmf(k = 1:3, project = "force")

  # The $id in admix_obj is the same as in the gen_tibble
  admix_obj$id

  # Reorder the q matrices based on the grouping variable
  admix_obj <- gt_admix_reorder_q(admix_obj,
    group = example_gt$population
  )

  # The $id in admix_obj is now reordered according to the population
  admix_obj$id
}

Run ADMIXTURE from R

Description

This function runs ADMIXTURE, taking either a gen_tibble or a file as an input. This is a wrapper that runs ADMIXTURE from the command line, and reads the output into R. It can run multiple values of k and multiple repeats for each k.

Usage

gt_admixture(
  x,
  k,
  n_runs = 1,
  crossval = FALSE,
  n_cores = 1,
  seed = NULL,
  conda_env = "auto",
  outdir = NULL
)

Arguments

Details

This is a wrapper for the command line program ADMIXTURE. It can either use a binary present in the main environment, or use a copy installed in a conda environment.

Value

an object of class gt_admix consisting of a list with the following elements:

• k the number of clusters

• Q a matrix with the admixture proportions

• P a matrix with the allele frequencies, with one row per locus in the same order as the loci of the input (as given by show_loci() for a gen_tibble, or the bim file for a bed file). ADMIXTURE reports the frequency of the second allele in the bim file, which for a gen_tibble is allele_ref.

• log a log of the output generated by ADMIXTURE (usually printed on the screen when running from the command line)

• cv the cross validation error (if crossval is TRUE)

• loglik the log likelihood of the model

• id the id column of the input gen_tibble (if applicable)

• group the group column of the input gen_tibble (if applicable)

References

Alexander, D.H., Novembre, J. and Lange, K. (2009) ‘Fast model-based estimation of ancestry in unrelated individuals’, Genome Research, 19(9), pp. 1655–1664. Available at: https://doi.org/10.1101/gr.094052.109.

Examples

# run the example only if we have the package installed
## Not run: 
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)
lobsters <- lobsters %>% group_by(population)
gt_admixture(lobsters,
  k = 2:3, seed = c(1, 2),
  n_runs = 2, crossval = TRUE
)

## End(Not run)

Convert a gen_tibble to a genind object from adegenet

Description

This function converts a gen_tibble to a genind object from adegenet. genind objects support arbitrary ploidy (including mixed-ploidy data across individuals), so this works on diploid, polyploid and mixed-ploidy gen_tibble objects alike.

Usage

gt_as_genind(x)

Arguments

Value

a genind object

Note

pseudohaploid individuals are treated as diploid (ploidy 2) using their stored genotype counts (dosage 0/2) directly, the same convention used by other tools (e.g. PLINK) when handling pseudohaploid data.

Examples

example_gt <- load_example_gt("gen_tbl")

# Convert to genind
gt_genind <- example_gt %>% gt_as_genind()

# Check object class
class(gt_genind)

Convert a gen_tibble to a genlight object from adegenet

Description

This function converts a gen_tibble to a genlight object from adegenet

Usage

gt_as_genlight(x)

Arguments

Details

The chromosome and position of each locus are stored in the genlight object, and the alleles are stored as "ref/alt" in loc.all (genlight counts the second allele, i.e. the alternate allele). As genlight only supports single letter alleles, loc.all is not set (with a warning) if any locus has a missing or multi-character allele.

Value

a genlight object

Examples

example_gt <- load_example_gt("gen_tbl")

# Convert to genlight
gt_genlight <- example_gt %>% gt_as_genlight()

# Check object class
class(gt_genlight)

Convert a gentibble to a .geno file for sNMF from the LEA package

Description

This function writes a .geno file from a gen_tibble. Unless a file path is given, a file with suffix .geno is written in the same location as the .rds and .bk files that underpin the gen_tibble.

Usage

gt_as_geno_lea(x, file = NULL)

Arguments

Details

The genotypes are read and written directly from the file backed matrix (one locus at a time), so this function also works on datasets that do not fit in memory.

Value

the path of the .geno file

Note

the .geno format only supports haploid or diploid data, so this function only works on diploid (or pseudohaploid) gen_tibble objects. Pseudohaploid data is simply treated as regular diploid genotype counts (dosage 0/2), the same convention used by other tools (e.g. PLINK).

See Also

LEA::geno()

Examples

example_gt <- load_example_gt("gen_tbl")

# Write a geno file
gt_as_geno_lea(example_gt, file = paste0(tempfile(), "_example.geno"))

Convert a gen_tibble to a data.frame compatible with hierfstat

Description

This function converts a gen_tibble to a data.frame formatted to be used by hierfstat functions.

Usage

gt_as_hierfstat(x)

Arguments

Value

a data.frame with a column 'pop' and further column representing the genotypes (with alleles recoded as 1 and 2)

Note

hierfstat only supports haploid or diploid data, so this function only works on diploid (or pseudohaploid) gen_tibble objects. Pseudohaploid data is simply treated as regular diploid genotype counts (dosage 0/2), the same convention used by other tools (e.g. PLINK).

Examples

example_gt <- load_example_gt("gen_tbl")

# Convert to hierfstat format
gt_hierfstat <- example_gt %>% gt_as_hierfstat()

# Check object class
class(gt_hierfstat)

Export a gen_tibble object to PLINK bed format

Description

This function exports all the information of a gen_tibble object into a PLINK bed, ped or raw file (and associated files, i.e. .bim and .fam for .bed; .fam for .ped).

Usage

gt_as_plink(
  x,
  file = NULL,
  type = c("bed", "ped", "raw"),
  overwrite = TRUE,
  chromosomes_as_int = FALSE
)

Arguments

Details

If the gen_tibble has been read in from vcf format, family.ID in the resulting plink files will be the same as sample.ID. If the gen_tibble has a grouping variable, this will be used as the family.ID in the resulting plink files (with several grouping variables, their values are pasted with ".", e.g. pop1.F); otherwise, the population is used as the family.ID. Parents, sex and phenotype are taken from the columns paternal_ID, maternal_ID, sex and phenotype, if present, as created by gen_tibble() when reading PLINK files (bed, ped or bigSNP rds): sex is written as 1 for "male", 2 for "female" and 0 (unknown) otherwise; phenotype as 1 for "control" and 2 for "case", and other values (e.g. a quantitative phenotype) as they are; missing parents are written as 0, and missing phenotypes as -9. Genotypes are read and written directly from the file backed matrix, so all formats can be written for datasets that do not fit in memory (but note that ped and raw files are much larger, and thus slower to write, than bed files).

Value

the path of the saved file

Note

PLINK's bed, ped and raw formats only support diploid data, so this function only works on diploid (or pseudohaploid) gen_tibble objects. Pseudohaploid data is simply treated as regular diploid genotype counts (dosage 0/2), the same convention used by PLINK itself.

Examples

example_gt <- load_example_gt("gen_tbl")

# Write a bed file
example_gt %>% gt_as_plink(type = "bed", file = paste0(tempfile(), "_plink"))

# Write a ped file
example_gt %>% gt_as_plink(type = "ped", file = paste0(tempfile(), "_plink"))

# Write a raw file
example_gt %>% gt_as_plink(type = "raw", file = paste0(tempfile(), "_plink"))

Convert a gen_tibble to a VCF

Description

This function write a VCF from a gen_tibble. The genotypes are read and written in chunks of loci, so that this works on datasets that do not fit in memory.

Usage

gt_as_vcf(x, file = NULL, chunk_size = NULL, overwrite = FALSE)

Arguments

Value

the path of the .vcf file

Examples

example_gt <- load_example_gt("gen_tbl")

# Write a vcf file
example_gt %>% gt_as_vcf()

Run K-clustering on principal components

Description

This function implements the clustering procedure used in Discriminant Analysis of Principal Components (DAPC, Jombart et al. 2010). This procedure consists in running successive K-means with an increasing number of clusters (k), after transforming data using a principal component analysis (PCA). For each model, several statistical measures of goodness of fit are computed, which allows to choose the optimal k using the function gt_cluster_pca_best_k(). See details for a description of how to select the optimal k and vignette("adegenet-dapc") for a tutorial.

Usage

gt_cluster_pca(
  x = NULL,
  n_pca = NULL,
  k_clusters = c(1, round(nrow(x$u)/10)),
  method = c("kmeans", "ward"),
  n_iter = 1e+05,
  n_start = 10,
  quiet = FALSE
)

Arguments

Value

a gt_cluster_pca object, which is a subclass of gt_pca with an additional element 'cluster', a list with elements:

• 'method' the clustering method (either kmeans or ward)

• 'n_pca' number of principal components used for clustering

• 'k' the k values explored by the function

• 'WSS' within sum of squares for each k

• 'AIC' the AIC for each k

• 'BIC' the BIC for each k

• 'groups' a list, with each element giving the group assignments for a given k

References

Jombart T, Devillard S and Balloux F (2010) Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genetics 11:94. doi:10.1186/1471-2156-11-94

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA object
pca <- gt_pca_partialSVD(lobsters)

# Run clustering on the first 10 PCs
gt_cluster_pca(
  x = pca,
  n_pca = 10,
  k_clusters = c(1, 5),
  method = "kmeans",
  n_iter = 1e5,
  n_start = 10,
  quiet = FALSE
)

# Alternatively, use method "ward"
gt_cluster_pca(
  x = pca,
  n_pca = 10,
  k_clusters = c(1, 5),
  method = "ward",
  quiet = FALSE
)

Find the best number of clusters based on principal components

Description

This function selects the best k value based on a chosen metric and criterion. It is equivalent to plotting the metric against the k values, and selecting the k that fulfills a given criterion (see details for an explanation of each criterion). This function simply adds an element 'best_k' to the gt_cluster_pca returned by gt_cluster_pca(). The choice can be over-ridden simply by assigning a different value to that element (e.g. for an object x and a desired k of 8, simply use x$best_k <- 8)

Usage

gt_cluster_pca_best_k(
  x,
  stat = c("BIC", "AIC", "WSS"),
  criterion = c("diffNgroup", "min", "goesup", "smoothNgoesup", "goodfit"),
  quiet = FALSE
)

Arguments

Details

The analysis of data simulated under various population genetics models (see reference) suggested an ad-hoc rule for the selection of the optimal number of clusters. First important result is that BIC seems more efficient than AIC and WSS to select the appropriate number of clusters (see example). The rule of thumb consists in increasing K until it no longer leads to an appreciable improvement of fit (i.e., to a decrease of BIC). In the most simple models (island models), BIC decreases until it reaches the optimal K, and then increases. In these cases, the best rule amounts to choosing the lowest K. In other models such as stepping stones, the decrease of BIC often continues after the optimal K, but is much less steep, so a change in slope can be taken as an indication of where the best k lies.

This function provides a programmatic way to select k. Note that it is highly recommended to look at the graph of BIC versus the numbers of clusters, to understand and validate the programmatic selection. The criteria available in this function are:

• "diffNgroup": differences between successive values of the summary statistics (by default, BIC) are split into two groups using a Ward's clustering method (see ?hclust), to differentiate sharp decrease from mild decreases or increases. The retained K is the one before the first group switch. This criterion appears to work well for island/hierarchical models, and decently for isolation by distance models, albeit with some instability. It can be confounded by an initial, very sharp decrease of the test statistics. IF UNSURE ABOUT THE CRITERION TO USE, USE THIS ONE.

• "min": the model with the minimum summary statistics (as specified by stat argument, BIC by default) is retained. Is likely to work for simple island model, using BIC. It is likely to fail in models relating to stepping stones, where the BIC always decreases (albeit by a small amount) as K increases. In general, this approach tends to over-estimate the number of clusters.

• "goesup": the selected model is the K after which increasing the number of clusters leads to increasing the summary statistics. Suffers from inaccuracy, since i) a steep decrease might follow a small 'bump' of increase of the statistics, and ii) increase might never happen, or happen after negligible decreases. Is likely to work only for clear-cut island models.

• "smoothNgoesup": a variant of "goesup", in which the summary statistics is first smoothed using a lowess approach. Is meant to be more accurate than "goesup" as it is less prone to stopping to small 'bumps' in the decrease of the statistics.

• "goodfit": another criterion seeking a good fit with a minimum number of clusters. This approach does not rely on differences between successive statistics, but on absolute fit. It selects the model with the smallest K so that the overall fit is above a given threshold.

Value

a 'gt_cluster_pca' object with an added element 'best_k'

References

Jombart T, Devillard S and Balloux F (2010) Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genetics 11:94. doi:10.1186/1471-2156-11-94

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA object
pca <- gt_pca_partialSVD(lobsters)

# Run clustering on the first 10 PCs
cluster_pca <- gt_cluster_pca(
  x = pca,
  n_pca = 10,
  k_clusters = c(1, 5),
  method = "kmeans",
  n_iter = 1e5,
  n_start = 10,
  quiet = FALSE
)

# Find best K through minimum BIC
cluster_pca <- gt_cluster_pca_best_k(cluster_pca,
  stat = "BIC",
  criterion = "min",
  quiet = FALSE
)
# Best K is stored in the object
cluster_pca$best_k

Discriminant Analysis of Principal Components for gen_tibble

Description

This function implements the Discriminant Analysis of Principal Components (DAPC, Jombart et al. 2010). This method describes the diversity between pre-defined groups. When groups are unknown, use gt_cluster_pca() to infer genetic clusters. See 'details' section for a succinct description of the method, and the vignette in the package adegenet ("adegenet-dapc") for a tutorial.

Usage

gt_dapc(x, pop = NULL, n_pca = NULL, n_da = NULL, loadings_by_locus = TRUE)

Arguments

Details

The Discriminant Analysis of Principal Components (DAPC) is designed to investigate the genetic structure of biological populations. This multivariate method consists in a two-steps procedure. First, genetic data are transformed (centred, possibly scaled) and submitted to a Principal Component Analysis (PCA). Second, principal components of PCA are submitted to a Linear Discriminant Analysis (LDA). A trivial matrix operation allows to express discriminant functions as linear combination of alleles, therefore allowing one to compute allele contributions. More details about the computation of DAPC are to be found in the indicated reference.

Results can be visualised with autoplot.gt_dapc(), see the help of that method for the available plots. There are also gt_dapc_tidiers for manipulating the results. For the moment, this function returns objects of class adegenet::dapc which are compatible with methods from adegenet; graphical methods for DAPC are documented in adegenet::scatter.dapc (see ?scatter.dapc). This is likely to change in the future, so make sure you do not rely on the objects remaining compatible.

This function aligns with the guidelines proposed by Thia (2023) for the standardized application of DAPC to genotype data. Our default settings are designed to follow these recommendations, so that the number of principal components (n_pca) defaults to the smaller of k-1 and the number of available principal components (where k is the number of populations or clusters), and the number of discriminant functions (n_da) is set to the minimum of k-1 and n_pca. The user can override these defaults by specifying the n_pca and n_da arguments, but caution is advised when adjusting n_pca to avoid potential overfitting. We recommend users consult these guidelines and consider their individual dataset to ensure best practices.

Note that there is no current method to predict scores for individuals not included in the original analysis. This is because we currently do not have a mechanism to store the pca information in the object, and that is needed for prediction.

Value

an object of class adegenet::dapc

References

Jombart T, Devillard S and Balloux F (2010) Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genetics 11:94. doi:10.1186/1471-2156-11-94 Thia, J. A. (2023). Guidelines for standardizing the application of discriminant analysis of principal components to genotype data. Molecular Ecology Resources, 23, 523–538. https://doi.org/10.1111/1755-0998.13706

See Also

gt_cluster_pca() gt_cluster_pca_best_k() adegenet::dapc()

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA object
pca <- gt_pca_partialSVD(lobsters)

# Run DAPC on the `gt_pca` object, providing `pop` as factor
populations <- as.factor(lobsters$population)
gt_dapc(pca, n_pca = 6, n_da = 2, pop = populations)

# Run clustering on the first 10 PCs
cluster_pca <- gt_cluster_pca(
  x = pca,
  n_pca = 10,
  k_clusters = c(1, 5),
  method = "kmeans",
  n_iter = 1e5,
  n_start = 10,
  quiet = FALSE
)

# Find best k
cluster_pca <- gt_cluster_pca_best_k(cluster_pca,
  stat = "BIC",
  criterion = "min"
)

# Run DAPC on the `gt_cluster_pca` object
gt_dapc(cluster_pca, n_pca = 10, n_da = 2)

#  should be stored (TRUE, default) or not (FALSE). This information is
#  required to predict group membership of new individuals using predict, but
#  makes the object slightly bigger.

Compute and store blocked f2 statistics for ADMIXTOOLS 2

Description

Compute and store blocked f2 statistics for ADMIXTOOLS 2

Usage

gt_extract_f2(
  .x,
  outdir = NULL,
  blgsize = 0.05,
  maxmem = 8000,
  maxmiss = 0,
  minmaf = 0,
  maxmaf = 0.5,
  minac2 = FALSE,
  outpop = NULL,
  outpop_scale = TRUE,
  transitions = TRUE,
  transversions = TRUE,
  overwrite = FALSE,
  adjust_pseudohaploid = NULL,
  fst = TRUE,
  afprod = TRUE,
  poly_only = c("f2"),
  apply_corr = TRUE,
  n_cores = 1,
  quiet = FALSE
)

Arguments

Value

SNP metadata (invisibly)

References

Maier R, Patterson N (2024). admixtools: Inferring demographic history from genetic data. R package version 2.0.4, https://github.com/uqrmaie1/admixtools.

This function prepares data for various ADMIXTOOLS 2 functions from the package ADMIXTOOLS 2. It takes a gen_tibble, computes allele frequencies and blocked f2-statistics, and writes the results to outdir. It is equivalent to admixtools::extract_f2().

Examples

bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)
lobsters <- lobsters %>% group_by(population)
f2_path <- tempfile()
gt_extract_f2(lobsters, outdir = f2_path, quiet = TRUE)
admixtools::f2_from_precomp(f2_path, verbose = FALSE)

Convert a genlight object from adegenet to a gen_tibble

Description

This function converts a genlight object from the adegenet package to a gen_tibble object

Usage

gt_from_genlight(x, backingfile = NULL, ...)

Arguments

Details

• Currently supports diploid genlight objects only (all values in ⁠@ploidy⁠ must be 2).

• Requires non-missing slots: loc.names, n.loc, loc.all, chromosome, position, ploidy, ind.names. The pop slot is optional; if absent, the returned gen_tibble will omit the population column.

Value

A gen_tibble object

Examples

# Create a simple genlight object
x <- new("genlight",
  list(
    indiv1 = c(1, 1, 0, 1, 1, 0),
    indiv2 = c(2, 1, 1, 0, 0, 0)
  ),
  ploidy = c(2, 2),
  loc.names = paste0("locus", 1:6),
  chromosome = c("chr1", "chr1", "chr2", "chr2", "chr3", "chr3"),
  position = c(100, 200, 150, 250, 300, 400),
  loc.all = c("A/T", "C/G", "G/C", "A/T", "T/C", "G/A"),
  pop = c("pop1", "pop2")
)


file <- paste0(tempfile(), "gt_from_genlight")
# Convert to gen_tibble
new_gt <- gt_from_genlight(x, backingfile = file)

Get the names of files storing the genotypes of a gen_tibble

Description

A function to return the names of the files used to store data in a gen_tibble. Specifically, this returns the .rds file storing the big

Usage

gt_get_file_names(x)

Arguments

Value

a character vector with the names and paths of the two files

Examples

example_gt <- load_example_gt("gen_tbl")

# To retrieve the names of and paths to the .bk and .rds files use:
gt_get_file_names(example_gt)

Checks if a gen_tibble has been imputed

Description

This function checks if a dataset has been imputed. Note that having imputation does not mean that the imputed values are used.

Usage

gt_has_imputed(x)

Arguments

Value

boolean TRUE or FALSE depending on whether the dataset has been imputed

Examples

example_gt <- load_example_gt("gen_tbl")

# The initial gen_tibble contains no imputed values
example_gt %>% gt_has_imputed()

# Now impute the gen_tibble
example_gt <- example_gt %>% gt_impute_simple()

# And we can check it has been imputed
example_gt %>% gt_has_imputed()

Simple imputation based on allele frequencies

Description

This function provides a very simple imputation algorithm for gen_tibble objects by using the mode, mean or sampling from the allele frequencies. Each locus is imputed independently (and thus linkage information is ignored).

Usage

gt_impute_simple(x, method = c("mode", "mean0", "random"), n_cores = 1)

Arguments

Details

This function is a wrapper around bigsnpr::snp_fastImputeSimple().

Value

a gen_tibble with imputed genotypes

See Also

bigsnpr::snp_fastImputeSimple() which this function wraps.

Examples

example_gt <- load_example_gt("gen_tbl")

# Impute the gen_tibble
example_gt <- example_gt %>% gt_impute_simple()

# And we can check it has been imputed
example_gt %>% gt_has_imputed()

Imputation based XGBoost

Description

This function provides a simple imputation algorithm for gen_tibble objects based on local XGBoost models.

Usage

gt_impute_xgboost(
  x,
  alpha = 1e-04,
  size = 200,
  p_train = 0.8,
  n_cor = nrow(x),
  seed = NA,
  n_cores = 1,
  append_error = TRUE
)

Arguments

Details

This function is a wrapper around bigsnpr::snp_fastImpute(). The error rates from the xgboost, if appended, can be retrieved with attr(x$genotypes, "imputed_errors") where x is the gen_tibble.

Value

a gen_tibble with imputed genotypes

See Also

bigsnpr::snp_fastImpute() which this function wraps.

Examples

example_gt <- load_example_gt("gen_tbl")

# Impute the gen_tibble
example_gt <- example_gt %>% gt_impute_xgboost()

# And we can check it has been imputed
example_gt %>% gt_has_imputed()

Load a gen_tibble

Description

Load a gen_tibble previously saved with gt_save(). If the .rds and .bk files have not been moved, they should be found automatically. If they were moved, use reattach_to to point to the .rds file (the .bk file needs to be in the same directory as the .rds file).

Usage

gt_load(file = NULL, reattach_to = NULL)

Arguments

Value

a gen_tibble

See Also

gt_save()

Examples

example_gt <- load_example_gt("gen_tbl")

# remove some individuals
example_gt_filtered <- example_gt %>% filter(id != "a")

# save the filtered gen_tibble object
backing_files <- gt_save(example_gt_filtered,
  file_name = paste0(tempfile(), "_example_filtered")
)

# backing_files[1] contains the name of the saved .gt file
backing_files[1]

# To load the saved gen_tibble object, use the path to the saved .gt file
reloaded_gt <- gt_load(backing_files[1])

# And we have loaded the gt without individual "a"
reloaded_gt

Order the loci table of a gen_tibble

Description

This function reorders the loci table so that positions within a chromosome are sequential. It also re-saves the genotypes into a new file backed matrix with the new order, so that it can be used by functions such as loci_ld_clump() and gt_pca_autoSVD(). If the loci table is already ordered, the original gen_tibble is returned. This function will update the backingfiles of the gen_tibble and return the gen_tibble object, use ⁠<-⁠ as per the example provided to ensure that the names of the newly updated backingfiles are stored in the gen_tibble object.

Usage

gt_order_loci(
  .x,
  use_current_table = FALSE,
  ignore_genetic_dist = TRUE,
  quiet = FALSE,
  ...
)

Arguments

Value

A gen_tibble

Examples

example_gt <- load_example_gt("gen_tbl") %>% select_loci(c(1, 5, 2, 6, 4, 3))

# Loci are in the wrong order
show_loci(example_gt)

# Reorder the loci, ignoring genetic distance
example_gt_ordered <- gt_order_loci(example_gt, ignore_genetic_dist = TRUE)

# Loci are now in the correct order
show_loci(example_gt_ordered)

Principal Component Analysis for gen_tibble objects

Description

There are a number of PCA methods available for gen_tibble objects. They are mostly designed to work on very large datasets, so they only compute a limited number of components. For smaller datasets, gt_pca_partialSVD() allows the use of partial (truncated) SVD to fit the PCA; this method is suitable when the number of individuals is much smaller than the number of loci. For larger dataset, gt_pca_randomSVD() is more appropriate. Finally, there is a method specifically designed for dealing with LD in large datasets, gt_pca_autoSVD(). Whilst this is arguably the best option, it is somewhat data hungry, and so only suitable for very large datasets (hundreds of individuals with several hundred thousands markers, or larger).

Note

Monomorphic markers must be removed before PCA is computed. The error message 'Error: some variables have zero scaling; remove them before attempting to scale.' indicates that monomorphic markers are present.

Using gt_pca_autoSVD with a small dataset will likely cause an error, see man page for details.

Rather than accessing elements of a gt_pca object directly, it is better to use tidy and augment. See gt_pca_tidiers.

See Also

Other gt_pca_functions: gt_pca_autoSVD(), gt_pca_partialSVD(), gt_pca_randomSVD()

PCA controlling for LD for gen_tibble objects

Description

This function performs Principal Component Analysis on a gen_tibble, using a fast truncated SVD with initial pruning and then iterative removal of long-range LD regions. This function is a wrapper for bigsnpr::snp_autoSVD()

Usage

gt_pca_autoSVD(
  x,
  k = 10,
  fun_scaling = bigsnpr::snp_scaleBinom(),
  thr_r2 = 0.2,
  use_positions = TRUE,
  size = 100/thr_r2,
  roll_size = 50,
  int_min_size = 20,
  alpha_tukey = 0.05,
  min_mac = 10,
  max_iter = 5,
  n_cores = 1,
  verbose = TRUE,
  total_var = TRUE
)

Arguments

Details

Using gt_pca_autoSVD requires a reasonably large dataset, as the function iteratively removes regions of long range LD. If you encounter: 'Error in rollmean(): Parameter 'size' is too large.', roll_size exceeds the number of variants on at least one of your chromosomes. Try reducing 'roll_size' to avoid this error.

Value

a gt_pca object, which is a subclass of bigSVD; this is an S3 list with elements: A named list (an S3 class "big_SVD") of

• d, the eigenvalues (singular values, i.e. as variances),

• u, the scores for each sample on each component (the left singular vectors)

• v, the loadings (the right singular vectors)

• center, the centering vector,

• scale, the scaling vector,

• method, a string defining the method (in this case 'autoSVD'),

• call, the call that generated the object.

• loci, the loci used after long range LD removal.

Note

rather than accessing these elements directly, it is better to use tidy and augment. See gt_pca_tidiers.

Monomorphic markers must be removed before PCA is computed. The error message 'Error: some variables have zero scaling; remove them before attempting to scale.' indicates that monomorphic markers are present.

See Also

bigsnpr::snp_autoSVD() which this function wraps.

Other gt_pca_functions: gt_pca, gt_pca_partialSVD(), gt_pca_randomSVD()

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

show_loci(lobsters)$chromosome <- "1"

# Create PCA object, including total variance
gt_pca_autoSVD(lobsters,
  k = 10,
  roll_size = 20,
  total_var = TRUE
)
# Change number of components and exclude total variance
gt_pca_autoSVD(lobsters,
  k = 5,
  roll_size = 20,
  total_var = FALSE
)

PCA for gen_tibble objects by partial SVD

Description

This function performs Principal Component Analysis on a gen_tibble, by partial SVD through the eigen decomposition of the covariance. It works well if the number of individuals is much smaller than the number of loci; otherwise, gt_pca_randomSVD() is a better option. This function is a wrapper for bigstatsr::big_SVD().

Usage

gt_pca_partialSVD(
  x,
  k = 10,
  fun_scaling = bigsnpr::snp_scaleBinom(),
  total_var = TRUE
)

Arguments

Details

With up to a few thousand individuals, this is usually much faster than gt_pca_randomSVD() (on 1,043 individuals and 643,733 loci, about 31 s on one core against about 175 s). The covariance between individuals is kept in memory, which takes 8 bytes per pair of individuals (about 0.2 GB for 5,000 individuals, and 3.2 GB for 20,000); with many more individuals, gt_pca_randomSVD() uses much less memory.

The genotypes are read in blocks of loci, decoded as doubles, of about 1 GB (set by options(bigstatsr.block.sizeGB)). A smaller size reduces the memory used, at the cost of some speed: on the same data, blocks of 0.25 GB (options(bigstatsr.block.sizeGB = 0.25)) reduced the peak memory from about 2.0 to 0.6 GB, and took about 25% longer.

Value

a gt_pca object, which is a subclass of bigSVD; this is an S3 list with elements: A named list (an S3 class "big_SVD") of

• d, the eigenvalues (singular values, i.e. as variances),

• u, the scores for each sample on each component (the left singular vectors)

• v, the loadings (the right singular vectors)

• center, the centering vector,

• scale, the scaling vector,

• method, a string defining the method (in this case 'partialSVD'),

• call, the call that generated the object.

• square_frobenius, used to compute the proportion of variance explained by the components (optional)

Note

Rather than accessing these elements directly, it is better to use tidy and augment. See gt_pca_tidiers.

Monomorphic markers must be removed before PCA is computed. The error message 'Error: some variables have zero scaling; remove them before attempting to scale.' indicates that monomorphic markers are present.

See Also

bigstatsr::big_SVD() which this function wraps.

Other gt_pca_functions: gt_pca, gt_pca_autoSVD(), gt_pca_randomSVD()

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA object, including total variance
gt_pca_partialSVD(lobsters,
  k = 10,
  total_var = TRUE
)
# Change number of components and exclude total variance
gt_pca_partialSVD(lobsters,
  k = 5,
  total_var = FALSE
)

PCA for gen_tibble objects by randomized partial SVD

Description

This function performs Principal Component Analysis on a gen_tibble, by randomised partial SVD based on the algorithm in RSpectra (by Yixuan Qiu and Jiali Mei).
This algorithm is linear in time in all dimensions and is very memory-efficient. Thus, it can be used on very large big.matrices. This function is a wrapper for bigstatsr::big_randomSVD()

Usage

gt_pca_randomSVD(
  x,
  k = 10,
  fun_scaling = bigsnpr::snp_scaleBinom(),
  tol = 1e-04,
  verbose = FALSE,
  n_cores = 1,
  fun_prod = bigstatsr::big_prodVec,
  fun_cprod = bigstatsr::big_cprodVec,
  total_var = TRUE
)

Arguments

Details

With up to a few thousand individuals, gt_pca_partialSVD() is usually much faster: on 1,043 individuals and 643,733 loci, this function took about 175 s on one core (with the default tol), against about 31 s for gt_pca_partialSVD(). This algorithm is most useful when there are too many individuals for gt_pca_partialSVD(), whose memory grows with the square of the number of individuals. Using several cores (n_cores), or a larger tol, also reduces its time.

Value

a gt_pca object, which is a subclass of bigSVD; this is an S3 list with elements: A named list (an S3 class "big_SVD") of

• d, the eigenvalues (singular values, i.e. as variances),

• u, the scores for each sample on each component (the left singular vectors)

• v, the loadings (the right singular vectors)

• center, the centering vector,

• scale, the scaling vector,

• method, a string defining the method (in this case 'randomSVD'),

• call, the call that generated the object.

Note

Rather than accessing these elements directly, it is better to use tidy and augment. See gt_pca_tidiers.

Monomorphic markers must be removed before PCA is computed. The error message 'Error: some variables have zero scaling; remove them before attempting to scale.' indicates that monomorphic markers are present.

See Also

bigstatsr::big_randomSVD() which this function wraps.

Other gt_pca_functions: gt_pca, gt_pca_autoSVD(), gt_pca_partialSVD()

Examples

vcf_path <-
  system.file("extdata", "anolis",
    "punctatus_t70_s10_n46_filtered.recode.vcf.gz",
    package = "tidypopgen"
  )
anole_gt <-
  gen_tibble(vcf_path, quiet = TRUE, backingfile = tempfile("anolis_"))

# Remove monomorphic loci and impute
anole_gt <- anole_gt %>% select_loci_if(loci_maf(genotypes) > 0)
anole_gt <- gt_impute_simple(anole_gt, method = "mode")

# Create PCA object, including total variance
gt_pca_randomSVD(anole_gt, k = 10, total_var = TRUE)

pcadapt analysis on a gen_tibble object

Description

pcadapt is an algorithm that detects genetic markers under selection. It is based on the principal component analysis (PCA) of the genotypes of the individuals. The method is described in Luu et al. (2017). See the R package pcadapt, which provides extensive documentation and examples.

Usage

gt_pcadapt(x, pca, k, n_cores = 1)

Arguments

Details

Internally, this function uses the snp_pcadapt function from the bigsnpr package.

Value

An object of subclass gt_pcadapt, a subclass of mhtest.

References

Luu, K., Bazin, E., Blum, M. G. B., & François, O. (2017). pcadapt: an R package for genome scans for selection based on principal component analysis. Molecular Ecology Resources, 17(1), 67–77.

See Also

bigsnpr::snp_pcadapt() which this function wraps.

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Create PCA object
pca <- gt_pca_partialSVD(lobsters)

# Create a gt_pcadapt object
gt_pcadapt(lobsters, pca, k = 2)

Set the ploidy of a gen_tibble to include pseudohaploids

Description

If a gen_tibble includes pseudohaploid data, its ploidy is set to -2 to indicate that some individuals are coded as pseudohaploids. The ploidy of the individuals is updated, with pseudohaploids set to 1 and diploids set to 2. However, the dosages are not changed, meaning that pseudohaploids are still coded as 0 or 2. If the gen_tibble is already set to pseudohaploid, running gt_pseudohaploid will update the ploidy values again, if pseudohaploid individuals have been removed then ploidy is reset to 2.

Usage

gt_pseudohaploid(x, test_n_loci = 10000)

Arguments

Value

a gen_tibble object with the ploidy set to -2 and the individual ploidy values updated to 1 or 2.

Examples

example_gt <- load_example_gt("gen_tbl")

# Detect pseudohaploids and set ploidy for the whole gen_tibble
example_gt <- example_gt %>% gt_pseudohaploid(test_n_loci = 3)

# Ploidy is now set to -2
show_ploidy(example_gt)

# Individual ploidy now varies between 1 (pseudohaploid) and 2 (diploid)
indiv_ploidy(example_gt)

Inbreeding from runs of homozygosity (F_ROH)

Description

Computes F_ROH, the proportion of the genome of each individual that is in runs of homozygosity, from the runs found by windows_indiv_roh(). F_ROH can be computed from the runs of at least a given length (e.g. 1, 2, 4, 8 and 16 Mb): long runs come from recent common ancestors, while short runs also reflect older inbreeding.

Usage

gt_roh_froh(
  roh,
  min_length_mb = 0,
  by = c("genome", "chromosome"),
  genome_length = NULL
)

Arguments

Details

The length of the genome is, by default, the sum of the lengths of the chromosomes, each being its largest position (as in detectRUNS); it can be given instead with genome_length (e.g. the length of the autosomes of the reference genome, or the length covered by the loci).

Value

a tibble with one row per individual (including those with no runs), with the grouping variables of roh (if any), id, the column min_length_mb if there are several values, chromosome if by = "chromosome", length_bp (the total length of the runs) and froh.

See Also

autoplot.indiv_roh() to plot F_ROH.

Examples

bed_file <- system.file("extdata", "roh", "roh_sim.bed",
  package = "tidypopgen"
)
sim_gt <- gen_tibble(bed_file, backingfile = tempfile(), quiet = TRUE)
sim_roh <- windows_indiv_roh(sim_gt %>% group_by(population))
# F_ROH of each individual
gt_roh_froh(sim_roh)
# from the runs of at least 1, 2 and 4 Mb
gt_roh_froh(sim_roh, min_length_mb = c(1, 2, 4))
# on each chromosome
gt_roh_froh(sim_roh, by = "chromosome")

Islands of runs of homozygosity

Description

Finds the regions of the genome where many individuals (at least a proportion threshold of them, overall or in each group) are in runs of homozygosity: islands of runs of homozygosity (also called ROH hotspots), which can point to regions under selection.

Usage

gt_roh_islands(roh, loci, threshold = 0.5, by_group = TRUE)

Arguments

Details

Islands are stretches of consecutive loci of a chromosome that are in runs in at least threshold of the individuals (see gt_roh_loci()).

Value

a tibble with one row per island, with the grouping variables (if counting by group), chromosome, start and end (the positions of its first and last loci), length_bp (end - start + 1), n_loci, first_locus and last_locus (rows of the loci table), and max_prop and mean_prop (the largest and mean proportions of individuals in a run over its loci)

Examples

bed_file <- system.file("extdata", "roh", "roh_sim.bed",
  package = "tidypopgen"
)
sim_gt <- gen_tibble(bed_file, backingfile = tempfile(), quiet = TRUE)
sim_roh <- windows_indiv_roh(sim_gt %>% group_by(population))
# regions in runs in at least 20% of the individuals of a population
gt_roh_islands(sim_roh, show_loci(sim_gt), threshold = 0.2)

The loci in runs of homozygosity

Description

Counts, for each locus, the individuals that are in a run of homozygosity (as found by windows_indiv_roh()), overall or in each group. Loci that are in runs in most individuals of a population (see gt_roh_islands()) can point to regions under selection. The counts are those of the .hom.summary file of PLINK --homozyg.

Usage

gt_roh_loci(roh, loci, by_group = TRUE)

Arguments

Value

a tibble with one row per locus (and per group), with the grouping variables (if counting by group), name, chromosome, position, n_roh (the number of individuals in a run), n_indiv (the number of individuals) and prop_roh (n_roh / n_indiv)

See Also

gt_roh_islands() to find the regions in runs in many individuals, and autoplot.indiv_roh() (type = "coverage") to plot them.

Examples

bed_file <- system.file("extdata", "roh", "roh_sim.bed",
  package = "tidypopgen"
)
sim_gt <- gen_tibble(bed_file, backingfile = tempfile(), quiet = TRUE)
sim_roh <- windows_indiv_roh(sim_gt %>% group_by(population))
roh_loci <- gt_roh_loci(sim_roh, show_loci(sim_gt))
# the loci most often in runs in each population
roh_loci %>%
  group_by(group) %>%
  slice_max(prop_roh, n = 1, with_ties = FALSE)

Detect runs of homozygosity with detectRUNS (deprecated)

Description

This function is deprecated: use windows_indiv_roh(), which finds runs of homozygosity as PLINK (--homozyg), with its parameters. It is kept so that old code still gives the same results: it calls the sliding-window algorithm of the package detectRUNS (detectRUNS::slidingRUNS.run()) on each individual, as windows_indiv_roh() used to do. It needs detectRUNS to be installed, and will be removed in a future version of tidypopgen.

Usage

gt_roh_window(
  .x,
  window_size = 15,
  threshold = 0.05,
  min_snp = 3,
  heterozygosity = FALSE,
  max_opp_window = 1,
  max_miss_window = 1,
  max_gap = 10^6,
  min_length_bps = 1000,
  min_density = 1/1000,
  max_opp_run = NULL,
  max_miss_run = NULL
)

Arguments

Details

If the gen_tibble is grouped, then the grouping variable is used to fill in the 'group' column (with several grouping variables, their values are pasted with ".", e.g. pop1.F). Otherwise, the 'group' column is filled with the same values as the 'id' column.

Value

A dataframe with RUNs of Homozygosity or Heterozygosity in the analysed dataset, in the format of detectRUNS, with columns "group", "id", "chrom", "nSNP", "from", "to", "lengthBps" (group: population, breed, case/control etc.; id: individual identifier; chrom: chromosome on which the run is located; nSNP: number of SNPs in the run; from: starting position of the run, in bps; to: end position of the run, in bps; lengthBps: size of the run)

References

Biscarini F, Cozzi P, Gaspa G, Marras G (2019) detectRUNS: Detect Runs of Homozygosity and Runs of Heterozygosity in Diploid Genomes. R package. doi:10.32614/CRAN.package.detectRUNS

See Also

detectRUNS::slidingRUNS.run() which this function wraps.

Examples

bed_file <- system.file("extdata", "roh", "roh_sim.bed",
  package = "tidypopgen"
)
sim_gt <- gen_tibble(bed_file, backingfile = tempfile(), quiet = TRUE)
# deprecated: use windows_indiv_roh() instead
sim_roh <- suppressWarnings(gt_roh_window(sim_gt %>% group_by(population)))
head(sim_roh)

Save a gen_tibble

Description

Save the tibble (and update the backing files). The gen_tibble object is saved to a file with extension .gt, together with update its .rds and .bk files. Note that multiple .gt files can be linked to the same .rds and .bk files; generally, this occurs when we create multiple subsets of the data. The .gt file then stores the information on what subset of the full dataset we are interested in, whilst the .rds and .bk file store the full dataset. To reload a gen_tibble, you can pass the name of the .gt file with gt_load().

Usage

gt_save(x, file_name = NULL, quiet = FALSE)

Arguments

Value

the file name and path of the .gt file, together with the .rds and .bk files

See Also

gt_load()

Examples

example_gt <- load_example_gt("gen_tbl")

# remove some individuals
example_gt <- example_gt %>% filter(id != "a")

# save filtered gen_tibble object
gt_save(example_gt, file_name = paste0(tempfile(), "_example_filtered"))

Sets a gen_tibble to use imputed data

Description

This function sets or unsets the use of imputed data. For some analysis, such as PCA, that does not allow for missing data, we have to use imputation, but for other analysis it might be preferable to allow for missing data.

Usage

gt_set_imputed(x, set = NULL)

Arguments

Value

the gen_tibble, invisibly

Examples

example_gt <- load_example_gt("gen_tbl")

# Impute the gen_tibble
example_gt <- example_gt %>% gt_impute_simple()

# Check whether the gen_tibble uses imputed values
example_gt %>% gt_uses_imputed()

# Set the gen_tibble to use imputed values
example_gt %>% gt_set_imputed(TRUE)

# And check that the gen_tibble uses imputed values again
example_gt %>% gt_uses_imputed()

Run SNMF from R in tidypopgen

Description

Run SNMF from R in tidypopgen

Usage

gt_snmf(
  x,
  k,
  project = "continue",
  n_runs = 1,
  alpha,
  tolerance = 1e-05,
  entropy = FALSE,
  percentage = 0.05,
  I,
  iterations = 200,
  ploidy = 2,
  seed = -1
)

Arguments

Details

This is a wrapper for LEA::snmf().

Value

an object of class gt_admix consisting of a list with the following elements:

• k the number of clusters

• Q a matrix with the admixture proportions

• P a matrix with the allele frequencies

• log a log of the output generated by ADMIXTURE (usually printed on the screen when running from the command line)

• cv the masked cross-entropy (if entropy is TRUE)

• loglik the log likelihood of the model

• id the id column of the input gen_tibble (if applicable)

• group the group column of the input gen_tibble (if applicable)

See Also

LEA::snmf()

Examples

# run the example only if we have the package installed
example_gt <- load_example_gt("gen_tbl")

# To run SNMF on a gen_tibble:
example_gt %>% gt_snmf(
  k = 1:3, project = "force", entropy = TRUE,
  percentage = 0.5, n_runs = 1, seed = 1, alpha = 100
)

Update the backing matrix

Description

This functions forces a re-write of the file backing matrix to match the gen_tibble. Individuals and loci are subsetted and reordered according to the current state of the gen_tibble. A .gt, .bk and .rds file will be created.

Usage

gt_update_backingfile(
  .x,
  backingfile = NULL,
  chunk_size = NULL,
  rm_unsorted_dist = TRUE,
  quiet = FALSE
)

Arguments

Details

This function does not check whether the positions of your genetic loci are sorted. To check this, and update the file backing matrix, use gt_order_loci(). Tests for this function are in test_gt_order_loci.R

Value

a gen_tibble with updated .gt, .bk, and .rds files (i.e. a new File Backed Matrix)

Examples

example_gt <- load_example_gt("gen_tbl")

# Here, new backingfiles are created, but without using `<-` the object
# loaded in the R session is not updated
gt_update_backingfile(example_gt)

# Make sure to use `<-` to update the file names associated with the
# gen_tibble object loaded in the R session
example_gt <- example_gt %>% gt_update_backingfile()

Checks if a gen_tibble uses imputed data

Description

This function checks if a dataset uses imputed data. Note that it is possible to have a dataset that has been imputed but it is currently not using imputation.

Usage

gt_uses_imputed(x)

Arguments

Value

boolean TRUE or FALSE depending on whether the dataset is using the imputed values

Examples

example_gt <- load_example_gt("gen_tbl")

# Impute the gen_tibble
example_gt <- example_gt %>% gt_impute_simple()

# Check whether the gen_tibble uses imputed values
example_gt %>% gt_uses_imputed()

Add dendrograms to a heatmap of pairwise distances

Description

This function adds dendrograms to a heatmap of pairwise distances, generated with heatmap_pairwise(). The heatmap must have been created with an order function based on hclust(), and the dendrograms will be based on the same clustering. Dendrograms can be added to the left (row) and/or top (column) of the heatmap. Note that the heatmap needs to be formatted first, as it can not be modified easily once the dendrogram panels have been added.

Usage

heatmap_add_dendro(
  plot,
  side = c("both", "left", "top"),
  rel_size = 0.15,
  line_color = "grey30",
  line_size = 0.5
)

Arguments

Details

To correctly align the dendrogram with the labels of the heatmap, we use a gtable to insert the dendrogram panel directly adjacent to the heatmap panel. This means that the dendrogram will be perfectly aligned with the heatmap cells regardless of axis label width. However, manipulating individual plots in a gtable is not easy, so the heatmap should be formatted fully before the dendrograms are added.

Value

A heatmap_dendro object (a subclass of gtable), which auto-draws on print() and is accepted by plot() and ggsave(), just like a ggplot or patchwork object.

Examples

example_gt <- load_example_gt("gen_tbl")
# With a pairwise matrix:
fst_dist <- example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Nei87", type = "pairwise")
mat_plot <- heatmap_pairwise(fst_dist,
  order = function(x) {
    stats::hclust(stats::as.dist(x))
  }
) +
  ggplot2::scale_fill_viridis_c()
heatmap_add_dendro(mat_plot)

Pairwise distance heatmap

Description

Plot a square heatmap from either a full symmetric distance matrix or a tidy tibble of distances. Rows and columns can optionally be reordered using either an explicit ordering vector or a function that computes an ordering (e.g. hierarchical clustering). Dendrograms can be added to the heatmap with heatmap_add_dendro(). Note that any decoration (e.g. axis labels, fill scale) should be added to the heatmap before adding dendrograms, as the latter will modify the plot structure and make it difficult to add further customization afterwards.

Usage

heatmap_pairwise(x, order = NULL)

Arguments

Value

A ggplot object.

See Also

heatmap_add_dendro()

Examples

example_gt <- load_example_gt("gen_tbl")

# With a pairwise matrix:
fst_dist <- example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Nei87", type = "pairwise")

heatmap_pairwise(fst_dist)

# Explicit ordering
heatmap_pairwise(
  fst_dist,
  order = rev(rownames(fst_dist))
) +
  ggplot2::scale_fill_viridis_c()

# Hierarchical clustering
heatmap_pairwise(
  fst_dist,
  order = function(x) {
    stats::hclust(stats::as.dist(x), method = "complete")
  }
)

# With a tidy tibble of distances:
fst_tidy <- example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Nei87", type = "tidy")

heatmap_pairwise(fst_tidy)

heatmap_pairwise(
  fst_tidy,
  order = function(x) {
    stats::hclust(stats::as.dist(x))
  }
)

Estimate individual observed heterozygosity

Description

Estimate observed heterozygosity (H_obs) for each individual (i.e. the frequency of loci that are heterozygous in an individual). A locus is heterozygous whenever the dosage of the alternate allele is neither zero nor the individual's own ploidy, so this works for diploid, polyploid, and mixed-ploidy gen_tibble objects alike (each individual is evaluated against its own ploidy, so individuals of differing ploidy can be mixed freely).

Usage

indiv_het_obs(.x, as_counts = FALSE, ...)

## S3 method for class 'tbl_df'
indiv_het_obs(.x, as_counts = FALSE, ...)

## S3 method for class 'vctrs_bigSNP'
indiv_het_obs(.x, as_counts = FALSE, ...)

Arguments

Value

either:

• a vector of heterozygosities, one per individuals in the gen_tibble

• a matrix with two columns, where the first is the number of heterozygous loci for each individual and the second is the number of missing values for each individual

Examples

example_gt <- load_example_gt("gen_tbl")

example_gt %>% indiv_het_obs()

# For observed heterozygosity as counts:
example_gt %>% indiv_het_obs(as_counts = TRUE)

Individual inbreeding coefficient

Description

This function calculates the inbreeding coefficient for each individual based on the beta estimate from Weir and Goudet (2017).

Usage

indiv_inbreeding(.x, method = c("WG17"), allele_sharing_mat = NULL, ...)

## S3 method for class 'tbl_df'
indiv_inbreeding(.x, method = c("WG17"), allele_sharing_mat = NULL, ...)

## S3 method for class 'vctrs_bigSNP'
indiv_inbreeding(.x, method = c("WG17"), allele_sharing_mat = NULL, ...)

## S3 method for class 'grouped_gen_tbl'
indiv_inbreeding(.x, method = c("WG17"), allele_sharing_mat = NULL, ...)

Arguments

Value

a numeric vector of inbreeding coefficients.

References

Weir, BS and Goudet J (2017) A Unified Characterization of Population Structure and Relatedness. Genetics (2017) 206:2085

Examples

example_gt <- load_example_gt("gen_tbl")

example_gt %>% indiv_inbreeding(method = "WG17")

Estimate individual missingness

Description

Estimate missingness for each individual (i.e. the frequency of missing genotypes in an individual).

Usage

indiv_missingness(.x, as_counts, block_size, ...)

## S3 method for class 'tbl_df'
indiv_missingness(
  .x,
  as_counts = FALSE,
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

## S3 method for class 'vctrs_bigSNP'
indiv_missingness(
  .x,
  as_counts = FALSE,
  block_size = bigstatsr::block_size(length(.x), 1),
  ...
)

Arguments

Value

a vector of missingness, one per individuals in the gen_tibble

Examples

example_gt <- load_example_gt("gen_tbl")

example_gt %>% indiv_missingness()

# For missingness as counts:
example_gt %>% indiv_missingness(as_counts = TRUE)

Return individual ploidy

Description

Returns the ploidy for each individual.

Usage

indiv_ploidy(.x, ...)

## S3 method for class 'tbl_df'
indiv_ploidy(.x, ...)

## S3 method for class 'vctrs_bigSNP'
indiv_ploidy(.x, ...)

Arguments

Value

a vector of ploidy, one per individuals in the gen_tibble

Examples

example_gt <- load_example_gt("gen_tbl")

example_gt %>% indiv_ploidy()

Test if the loci table is ordered

Description

This functions checks that all SNPs in a chromosome are adjacent in the loci table, and that positions (and, unless ignore_genetic_dist = TRUE, genetic distances) are sorted within chromosomes. Missing positions (or missing genetic distances, if they are checked) are an error, as the order of those loci is unknown.

Usage

is_loci_table_ordered(
  .x,
  error_on_false = FALSE,
  ignore_genetic_dist = TRUE,
  ...
)

## S3 method for class 'tbl_df'
is_loci_table_ordered(
  .x,
  error_on_false = FALSE,
  ignore_genetic_dist = TRUE,
  ...
)

## S3 method for class 'vctrs_bigSNP'
is_loci_table_ordered(
  .x,
  error_on_false = FALSE,
  ignore_genetic_dist = TRUE,
  ...
)

Arguments

Value

TRUE if the loci are ordered, FALSE otherwise

Examples

example_gt <- load_example_gt("gen_tbl")

example_gt %>% is_loci_table_ordered()

Compute the Population Branch Statistics with block jackknife standard errors

Description

The function computes the genome-wide population branch statistics (PBS) for each combination of populations, together with its standard error estimated with a weighted block jackknife. The PBS is a measure of the genetic differentiation between one focal population and two reference populations.

Usage

jackknife_nwise_pop_pbs(.x, ...)

## S3 method for class 'gen_tbl'
jackknife_nwise_pop_pbs(
  .x,
  type = c("matrix", "tidy"),
  fst_method = c("Hudson", "Nei87", "WC84"),
  return_fst = FALSE,
  block_size,
  size_unit = c("snp", "cM", "bp"),
  ...
)

Arguments

Details

The genome is divided into non-overlapping blocks, which never span more than one chromosome, and pairwise Fst (and thus PBS) is recomputed leaving out one block at a time. The standard error is computed with the weighted delete-m jackknife of Busing et al. (1999), where each block is weighted by the number of loci it contains. See jackknife_stats_generic() for details.

Value

either a data frame with the following columns:

• quantity: the quantity reported in the row, one of estimate (computed on all loci), jack_estimate (the bias corrected jackknife estimate), se (the jackknife standard error), z (jack_estimate/se), and n_blocks (the number of blocks used)

• fst_a.b: the Fst value for population a and b, if return_fst is TRUE

• pbs_a.b.c: the PBS value for population a given b & c (there will be multiple such columns covering all 3 way combinations of populations in the grouped gen_tibble object)

• pbsn1_a.b.c: the normalised PBS value for population a given b & c or a tidy tibble with the following columns:

• stat_name: the name of populations used in the pbs calculation (e.g. "pbs_pop1.pop2.pop3"). If return_fst is TRUE, stat_name will also include "fst" calculations in the same column (e.g. "fst_pop1.pop2").

• estimate, jack_estimate, se, z, n_blocks: as described above

References

Busing FMTA, Meijer E, Van Der Leeden R (1999) Delete-m jackknife for unequal m. Statistics and Computing 9:3-8.

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# Calculate nwise pbs with blocks of 2 SNPs
example_gt %>%
  jackknife_nwise_pop_pbs(
    block_size = 2,
    size_unit = "snp"
  )

Compute pairwise Fst with block jackknife standard errors

Description

This function computes genome-wide pairwise Fst, together with its standard error estimated with a weighted block jackknife.

Usage

jackknife_pairwise_pop_fst(.x, ...)

## S3 method for class 'gen_tbl'
jackknife_pairwise_pop_fst(
  .x,
  type = c("matrix", "tidy"),
  method = c("Hudson", "Nei87", "WC84"),
  block_size,
  size_unit = c("snp", "cM", "bp"),
  ...
)

Arguments

Details

The genome is divided into non-overlapping blocks, which never span more than one chromosome, and Fst is recomputed leaving out one block at a time. The standard error is computed with the weighted delete-m jackknife of Busing et al. (1999), where each block is weighted by the number of loci it contains. Fst is computed as the ratio of the sum of the per locus numerators and denominators, so leaving out a block is achieved by subtracting its sums, without having to rescan the genotypes. See jackknife_stats_generic() for details.

Value

either a data frame with the following columns:

• quantity: the quantity reported in the row, one of estimate (Fst computed on all loci), jack_estimate (the bias corrected jackknife estimate), se (the jackknife standard error), z (jack_estimate/se), and n_blocks (the number of blocks used)

• fst_a.b: the values for the pairwise Fst between population a and b (there will be multiple such columns if there are more than two populations) or a tidy tibble with the following columns:

• stat_name: the name of population a and b used in the pairwise Fst calculation (e.g. "fst_pop1.pop2")

• estimate, jack_estimate, se, z, n_blocks: as described above

References

Busing FMTA, Meijer E, Van Der Leeden R (1999) Delete-m jackknife for unequal m. Statistics and Computing 9:3-8.

Examples

example_gt <- load_example_gt("gen_tbl")

example_gt %>%
  group_by(population) %>%
  jackknife_pairwise_pop_fst(
    block_size = 2,
    size_unit = "snp"
  )

Compute Tajima's D with block jackknife standard errors

Description

This function computes genome-wide Tajima's D, together with its standard error estimated with a weighted block jackknife.

Usage

jackknife_pop_tajimas_d(.x, ...)

## S3 method for class 'gen_tbl'
jackknife_pop_tajimas_d(
  .x,
  type = c("matrix", "tidy", "list"),
  block_size,
  size_unit = c("snp", "cM", "bp"),
  ...
)

Arguments

Details

The genome is divided into non-overlapping blocks, which never span more than one chromosome, and Tajima's D is recomputed leaving out one block at a time. The standard error is computed with the weighted delete-m jackknife of Busing et al. (1999), where each block is weighted by the number of loci it contains. Tajima's D is a function of sums across loci (of the nucleotide diversity, the number of segregating sites, and quantities that depend on the sample size of each locus), so leaving out a block is achieved by subtracting its sums, without having to rescan the genotypes. See jackknife_stats_generic() for details.

Missing genotypes are accounted for by using the number of non-missing allele copies at each locus as its sample size, following Bailey et al. (2025) as implemented in pixy. See pop_tajimas_d() for details. Loci with no genotyped individuals are ignored. The variance of Tajima (1989) is computed for the mean sample size across all loci (rounded to the nearest integer), and the same sample size is used when leaving out each block: recomputing the rounded mean for each block left out would make Tajima's D jump whenever this changes the rounding (e.g. when the mean sample size is close to x.5), which the jackknife would amplify into an inflated standard error and a biased jackknife estimate.

Note that the variance of Tajima (1989) assumes that there is no recombination among loci. Over a whole genome, Tajima's D is thus not standardised: it is approximately proportional to \hat{\theta}_\pi / \hat{\theta}_W - 1, with a factor that depends on the sample size (about 5 for 10 allele copies, 4.2 for 20, and 3.4 for 100), and with genome-wide data its standard error is much smaller than 1. Use the jackknife standard error, rather than the rules of thumb used for single loci, to judge whether genome-wide Tajima's D differs from 0, and compare groups with different sample sizes with caution. The z score tests whether genome-wide Tajima's D differs from 0, which most demographic histories other than a constant population size will cause. Blocks should be larger than the extent of linkage disequilibrium, or the standard error will be underestimated.

Value

if data is not grouped, a data frame with the following columns:

• estimate: Tajima's D computed on all loci

• jack_estimate: the bias corrected jackknife estimate

• se: the jackknife standard error

• z: the z score (jack_estimate/se)

• n_blocks: the number of blocks used if data are grouped, either: a data frame with the following columns:

• quantity: the quantity reported in the row (estimate, jack_estimate, se, z or n_blocks, as described above)

• group: the value of the quantity for the group (there will be as many of these columns as groups in the gen_tibble, and they will be named by the grouping levels) a tidy tibble with the following columns:

• group: the name of the group

• estimate, jack_estimate, se, z, n_blocks: as described above or a list of data frames, one per group, with the same columns as for ungrouped data.

References

Busing FMTA, Meijer E, Van Der Leeden R (1999) Delete-m jackknife for unequal m. Statistics and Computing 9:3-8.

Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123:585-595.

Bailey N, Stevison L, Samuk K (2025) Correcting for bias in estimates of \theta_W and Tajima's D from missing data in next-generation sequencing. Molecular Ecology Resources, e14104. doi:10.1111/1755-0998.14104

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# Calculate Tajima's D using blocks of 2 SNPs
example_gt %>%
  jackknife_pop_tajimas_d(
    block_size = 2,
    size_unit = "snp"
  )

Estimate a block jackknife standard error from per locus estimates

Description

This function is mostly designed for developers: it is a general function to estimate the genome-wide value of a statistic, together with its standard error, using a weighted block jackknife. It takes a vector of per locus estimates, and aggregates them by their mean (or a custom function), leaving out one block of loci at a time. To compute specific quantities directly from a gen_tibble, use the appropriate ⁠jackknife_*⁠ functions, e.g jackknife_pairwise_pop_fst() to compute pairwise Fst.

Usage

jackknife_stats_generic(
  .x,
  loci_table,
  operator = c("mean", "custom"),
  block_size,
  size_unit = c("snp", "cM", "bp"),
  f = NULL,
  ...
)

Arguments

Details

The genome is divided into non-overlapping blocks, which never span more than one chromosome. Blocks can be defined as a number of SNPs, a number of base pairs, or a size in centiMorgans (which requires the 'genetic_dist' column of the loci table to be filled in). Blocks in base pairs cover positions 1 to block_size, block_size + 1 to 2 * block_size, and so on along each chromosome (as the windows of windows_stats_generic() with step_size equal to window_size), and blocks in centiMorgans start from a genetic distance of 0. The statistic is then recomputed leaving out one block at a time, and the standard error is computed with the weighted delete-m jackknife of Busing et al. (1999), where each block is weighted by the number of loci with non-missing values that it contains. Blocks should be larger than the extent of linkage disequilibrium (e.g. 5 cM, or a few Mb).

The jackknife assumes that the statistic estimates a quantity that does not change with the number of loci (such as a mean, or a ratio of sums), and that is a smooth function of the data. A custom function f should thus not be a sum (leaving out a block always makes a sum smaller, which the jackknife would mistake for a bias), nor a statistic such as the median or other quantiles, for which the jackknife does not give valid standard errors. For the sum of the per locus estimates, use operator = "mean", and multiply the estimates and the standard error by the number of non-missing values.

Value

A tibble with columns: 'estimate' (the statistic computed on all loci), 'jack_estimate' (the bias corrected jackknife estimate), 'se' (the jackknife standard error), 'z' (the z score, jack_estimate/se), and 'n_blocks' (the number of blocks used).

References

Busing FMTA, Meijer E, Van Der Leeden R (1999) Delete-m jackknife for unequal m. Statistics and Computing 9:3-8.

Examples

example_gt <- load_example_gt("gen_tbl")

miss_by_locus <- loci_missingness(example_gt)

# Calculate the mean missingness and its standard error,
# using blocks of 2 SNPs
jackknife_stats_generic(miss_by_locus,
  loci_table = show_loci(example_gt),
  operator = "mean", block_size = 2,
  size_unit = "snp"
)

Load example gen_tibble

Description

This function creates a gen_tibble object for use in examples in documentation.

Usage

load_example_gt(
  type = c("gen_tbl", "grouped_gen_tbl", "grouped_gen_tbl_sf", "gen_tbl_sf",
    "gen_tbl_gen_dist")
)

Arguments

Value

an example object of the class gen_tbl.

Examples

# This function creates an example gen_tibble object
example_gt <- load_example_gt("gen_tbl")

Estimate allele frequencies at each locus

Description

Allele frequencies can be estimates as minimum allele frequencies (MAF) with loci_maf() or the frequency of the alternate allele (with loci_alt_freq()). The latter are in line with the genotypes matrix (e.g. as extracted by show_loci()). Most users will be in interested in the MAF, but the raw frequencies might be useful when computing aggregated statistics. Both loci_maf() and loci_alt_freq() have efficient methods to support grouped gen_tibble objects. These can return a tidied tibble, a list, or a matrix.

Usage

loci_alt_freq(
  .x,
  .col = "genotypes",
  as_counts = FALSE,
  n_cores,
  block_size,
  type,
  ...
)

## S3 method for class 'tbl_df'
loci_alt_freq(
  .x,
  .col = "genotypes",
  as_counts = FALSE,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

## S3 method for class 'vctrs_bigSNP'
loci_alt_freq(
  .x,
  .col = "genotypes",
  as_counts = FALSE,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(length(.x), 1),
  ...
)

## S3 method for class 'grouped_gen_tbl'
loci_alt_freq(
  .x,
  .col = "genotypes",
  as_counts = FALSE,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  type = c("tidy", "list", "matrix"),
  ...
)

loci_maf(.x, .col = "genotypes", n_cores, block_size, type, ...)

## S3 method for class 'tbl_df'
loci_maf(
  .x,
  .col = "genotypes",
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

## S3 method for class 'vctrs_bigSNP'
loci_maf(
  .x,
  .col = "genotypes",
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(length(.x), 1),
  ...
)

## S3 method for class 'grouped_gen_tbl'
loci_maf(
  .x,
  .col = "genotypes",
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  type = c("tidy", "list", "matrix"),
  ...
)

Arguments

Value

a vector of frequencies, one per locus, if as_counts = FALSE; else a matrix of two columns, the count of alternate alleles and the count valid alleles (i.e. the sum of alternate and reference)

Examples

example_gt <- load_example_gt("gen_tbl")

# For alternate allele frequency
example_gt %>% loci_alt_freq()

# For alternate allele frequency per locus per population
example_gt %>%
  group_by(population) %>%
  loci_alt_freq()
# alternatively, return a list of populations with their frequencies
example_gt %>%
  group_by(population) %>%
  loci_alt_freq(type = "list")
# or a matrix with populations in columns and loci in rows
example_gt %>%
  group_by(population) %>%
  loci_alt_freq(type = "matrix")
# or within reframe (not recommended, as it much less efficient
# than using it directly as shown above)
library(dplyr)
example_gt %>%
  group_by(population) %>%
  reframe(alt_freq = loci_alt_freq(genotypes))
# For MAF
example_gt %>% loci_maf()

# For minor allele frequency per locus per population
example_gt %>%
  group_by(population) %>%
  loci_maf()
# alternatively, return a list of populations with their frequencies
example_gt %>%
  group_by(population) %>%
  loci_maf(type = "list")
# or a matrix with populations in columns and loci in rows
example_gt %>%
  group_by(population) %>%
  loci_maf(type = "matrix")

Get the chromosomes of loci in a gen_tibble

Description

Extract the loci chromosomes from a gen_tibble (or directly from its genotype column).

Usage

loci_chromosomes(.x, .col = "genotypes", ...)

## S3 method for class 'tbl_df'
loci_chromosomes(.x, .col = "genotypes", ...)

## S3 method for class 'vctrs_bigSNP'
loci_chromosomes(.x, .col = "genotypes", ...)

Arguments

Value

a character vector of chromosomes

Examples

example_gt <- load_example_gt("gen_tbl")
example_gt %>% loci_chromosomes()

Test Hardy-Weinberg equilibrium at each locus

Description

Return the p-value from an exact test of HWE.

Usage

loci_hwe(.x, .col = "genotypes", ...)

## S3 method for class 'tbl_df'
loci_hwe(
  .x,
  .col = "genotypes",
  mid_p = TRUE,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

## S3 method for class 'vctrs_bigSNP'
loci_hwe(
  .x,
  .col = "genotypes",
  mid_p = TRUE,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(length(.x), 1),
  ...
)

## S3 method for class 'grouped_gen_tbl'
loci_hwe(
  .x,
  .col = "genotypes",
  mid_p = TRUE,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  type = c("tidy", "list", "matrix"),
  ...
)

Arguments

Details

This function uses the original C++ algorithm from PLINK 1.90.

Value

a vector of probabilities from HWE exact test, one per locus. For a grouped gen_tibble, empty groups (e.g. from group_by(.drop = FALSE)) are NA.

Author(s)

the C++ algorithm was written by Christopher Chang for PLINK 1.90, based on original code by Jan Wigginton (the code was released under GPL3).

Examples

example_gt <- load_example_gt("gen_tbl")

# For HWE
example_gt %>% loci_hwe()

# For loci_hwe per locus per population, use reframe
example_gt %>%
  group_by(population) %>%
  reframe(loci_hwe = loci_hwe(genotypes))

Clump loci based on a Linkage Disequilibrium threshold

Description

This function uses clumping to remove SNPs at high LD. When used with its default options, clumping based on MAF is similar to standard pruning (as done by PLINK with "–indep-pairwise (size+1) 1 thr.r2", but it results in a better spread of SNPs over the chromosome. This function is a wrapper around bigsnpr::snp_clumping(). See https://privefl.github.io/bigsnpr/articles/pruning-vs-clumping.html for more information on the differences between pruning and clumping.

Usage

loci_ld_clump(.x, .col = "genotypes", ...)

## S3 method for class 'tbl_df'
loci_ld_clump(.x, .col = "genotypes", ...)

## S3 method for class 'vctrs_bigSNP'
loci_ld_clump(
  .x,
  .col = "genotypes",
  S = NULL,
  thr_r2 = 0.2,
  size = 100/thr_r2,
  exclude = NULL,
  use_positions = TRUE,
  n_cores = 1,
  return_id = FALSE,
  ...
)

Arguments

Details

Any missing values in the genotypes of a gen_tibble passed to loci_ld_clump will cause an error. To deal with missingness, see gt_impute_simple().

Value

a boolean vector indicating whether the SNP should be kept (if 'return_id = FALSE', the default), else a vector of SNP indices to be kept (if 'return_id = TRUE')

See Also

bigsnpr::snp_clumping() which this function wraps.

Examples

example_gt <- load_example_gt("gen_tbl") %>% gt_impute_simple()

# To return a boolean vector indicating whether the SNP should be kept
example_gt %>% loci_ld_clump()
# To return a vector of SNP indices to be kept
example_gt %>% loci_ld_clump(return_id = TRUE)

Estimate missingness at each locus

Description

Estimate the rate of missingness at each locus. This function has an efficient method to support grouped gen_tibble objects, which can return a tidied tibble, a list, or a matrix.

Usage

loci_missingness(
  .x,
  .col = "genotypes",
  as_counts = FALSE,
  n_cores = bigstatsr::nb_cores(),
  block_size,
  type,
  ...
)

## S3 method for class 'tbl_df'
loci_missingness(
  .x,
  .col = "genotypes",
  as_counts = FALSE,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

## S3 method for class 'vctrs_bigSNP'
loci_missingness(
  .x,
  .col = "genotypes",
  as_counts = FALSE,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(length(.x), 1),
  ...
)

## S3 method for class 'grouped_gen_tbl'
loci_missingness(
  .x,
  .col = "genotypes",
  as_counts = FALSE,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  type = c("tidy", "list", "matrix"),
  ...
)

Arguments

Value

a vector of frequencies, one per locus

Examples

example_gt <- load_example_gt("gen_tbl")

# For missingness
example_gt %>% loci_missingness()

# For missingness per locus per population
example_gt %>%
  group_by(population) %>%
  loci_missingness()
# alternatively, return a list of populations with their missingness
example_gt %>%
  group_by(population) %>%
  loci_missingness(type = "list")
# or a matrix with populations in columns and loci in rows
example_gt %>%
  group_by(population) %>%
  loci_missingness(type = "matrix")
# or within reframe (not recommended, as it much less efficient
# than using it directly as shown above)
example_gt %>%
  group_by(population) %>%
  reframe(missing = loci_missingness(genotypes))

Get the names of loci in a gen_tibble

Description

Extract the loci names from a gen_tibble (or directly from its genotype column).

Usage

loci_names(.x, .col = "genotypes", ...)

## S3 method for class 'tbl_df'
loci_names(.x, .col = "genotypes", ...)

## S3 method for class 'vctrs_bigSNP'
loci_names(.x, .col = "genotypes", ...)

Arguments

Value

a character vector of names

Examples

example_gt <- load_example_gt("gen_tbl")
example_gt %>% loci_names()

Estimate nucleotide diversity (pi) at each locus

Description

Estimate nucleotide diversity (\pi) at each locus, accounting for missing values.

Usage

loci_pi(.x, .col = "genotypes", n_cores, block_size, type, ...)

## S3 method for class 'tbl_df'
loci_pi(
  .x,
  .col = "genotypes",
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

## S3 method for class 'vctrs_bigSNP'
loci_pi(
  .x,
  .col = "genotypes",
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(length(.x), 1),
  ...
)

## S3 method for class 'grouped_gen_tbl'
loci_pi(
  .x,
  .col = "genotypes",
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  type = c("tidy", "list", "matrix"),
  ...
)

Arguments

Details

Nucleotide diversity is estimated as

\pi = \frac{c_0 c_1}{n(n - 1)/2},

where:

• c_0 is the number of observed copies of allele 0,

• c_1 is the number of observed copies of allele 1,

• n = c_0 + c_1 is the total number of non-missing allele copies at the locus.

Missing values are accommodated by computing n, c_0, and c_1 using only non-missing genotype data. For polyploid data, each individual contributes a number of allele copies equal to its ploidy.

This corresponds to the unbiased gene diversity estimator of Nei and Roychoudhury (1974) and is applicable to diploid, polyploid mixed-ploidy, and pseudohaploid gen_tibble objects alike.

Value

a vector of frequencies, one per locus

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# For pi
example_gt %>% loci_pi()

# For pi per locus per population
example_gt %>%
  group_by(population) %>%
  loci_pi()
# alternatively, return a list of populations with their pi
example_gt %>%
  group_by(population) %>%
  loci_pi(type = "list")
# or a matrix with populations in columns and loci in rows
example_gt %>%
  group_by(population) %>%
  loci_pi(type = "matrix")
# or within reframe (not recommended, as it much less efficient
# than using it directly as shown above)
example_gt %>%
  group_by(population) %>%
  reframe(pi = loci_pi(genotypes))

Find transitions

Description

Use the loci table to define which loci are transitions

Usage

loci_transitions(.x, .col = "genotypes", ...)

## S3 method for class 'tbl_df'
loci_transitions(.x, .col = "genotypes", ...)

## S3 method for class 'vctrs_bigSNP'
loci_transitions(.x, .col = "genotypes", ...)

Arguments

Value

a logical vector defining which loci are transitions

Examples

example_gt <- load_example_gt("gen_tbl")
example_gt %>% loci_transitions()

Find transversions

Description

Use the loci table to define which loci are transversions

Usage

loci_transversions(.x, .col = "genotypes", ...)

## S3 method for class 'tbl_df'
loci_transversions(.x, .col = "genotypes", ...)

## S3 method for class 'vctrs_bigSNP'
loci_transversions(.x, .col = "genotypes", ...)

Arguments

Value

a logical vector defining which loci are transversions

Examples

example_gt <- load_example_gt("gen_tbl")
example_gt %>% loci_transversions()

Matrix with integer row/column names class

Description

Create a matrix with integer or character names for rows and columns. Integer names are useful when dealing with very large numbers of rows or columns, as the integer names can be stored more efficiently than character names. Whilst this is possible for data.frames, the native matrix class only allows for character row and column names via the dimnames attribute. Use row_names() and col_names() to get or set the integer names, which are stored in special attributes (int_rownames and int_colnames) on the matrix.

Usage

matrix_int_names(data, row_names = NULL, col_names = NULL)

Arguments

Details

This class allows you to have integer names as attributes while still being a matrix, and provides methods for getting and setting these names. When you create a matrix_int_names object, you can provide integer or character vectors for row and column names. If you provide integer vectors, they will be stored in special attributes (int_rownames and int_colnames) instead of the standard dimnames. If you provide character vectors, they will be stored in the standard dimnames as usual. You can also mix and match, having integer names for rows and character names for columns, or vice versa. Note that, since the row and column names are stored in special attributes, you have to use row_names() and col_names() to get or set them, rather than rownames() and colnames(), which will only return character names if present.

Value

A matrix_int_names object

See Also

Other matrix_int_names_functions: [.matrix_int_names(), as.data.frame.matrix_int_names(), col_names(), print.matrix_int_names(), row_names()

Examples

# Create a matrix with integer row and column names
my_mat <- matrix_int_names(matrix(1:6, nrow = 2),
  row_names = c(10L, 20L),
  col_names = c(100L, 200L, 300L)
)
row_names(my_mat) # returns integer row names
col_names(my_mat) # returns integer column names

A mutate method for gen_tibble objects

Description

A mutate method for gen_tibble objects

Usage

## S3 method for class 'gen_tbl'
mutate(..., deparse.level = 1)

Arguments

Value

a gen_tibble

Examples

example_gt <- load_example_gt("gen_tbl")

# Add a new column
example_gt %>% mutate(region = "East")

A mutate method for grouped gen_tibble objects

Description

A mutate method for grouped gen_tibble objects

Usage

## S3 method for class 'grouped_gen_tbl'
mutate(..., deparse.level = 1)

Arguments

Value

a grouped gen_tibble

Examples

test_gt <- load_example_gt("grouped_gen_tbl")
test_gt %>% mutate(region = "East")
test_gt <- load_example_gt("grouped_gen_tbl_sf")
test_gt %>% mutate(region = "East")

Compute the Population Branch Statistics for each combination of populations

Description

The function computes the population branch statistics (PBS) for each combination of populations at each locus. The PBS is a measure of the genetic differentiation between one focal population and two reference populations, and is used to identify outlier loci that may be under selection.

Usage

nwise_pop_pbs(.x, ...)

## S3 method for class 'gen_tbl'
nwise_pop_pbs(
  .x,
  type = c("tidy", "matrix"),
  fst_method = c("Hudson", "Nei87", "WC84"),
  return_fst = FALSE,
  ...
)

Arguments

Value

Either a matrix_int_names with locus ID as row names (use row_names() to retrieve them, as they can be integers) and the following columns:

• pbs_a.b.c: the PBS value for population a given b & c (there will be multiple such columns covering all 3 way combinations of populations in the grouped gen_tibble object)

• pbsn1_a.b.c: the normalized PBS value for population a given b & c.

• fst_a.b: the Fst value for population a and b, if return_fst is TRUE or a tidy tibble with the following columns:

• loci: the locus ID

• stat_name: the name of populations used in the pbs calculation (e.g. "pbs_pop1.pop2.pop3"). If return_fst is TRUE, stat_name will also include "fst" calculations in the same column (e.g. "fst_pop1.pop2").

• value: the pbs value for the populations

References

Yi X, et al. (2010) Sequencing of 50 human exomes reveals adaptation to high altitude. Science 329: 75-78.

Examples

example_gt <- load_example_gt()

# We can compute the PBS for all populations using "Hudson" method
example_gt %>%
  group_by(population) %>%
  nwise_pop_pbs(fst_method = "Hudson")

Compute the joint site frequency spectrum of all populations

Description

Compute the multidimensional (joint) site frequency spectrum (SFS) of all the populations of a grouped gen_tibble, i.e. the number of loci with each combination of numbers of copies of the derived (or minor) allele in the populations (the "multiSFS" of easySFS and fastsimcoal2). Missing genotypes are dealt with as in easySFS, by projecting each population down to n_proj allele copies (see pop_sfs() for details). Only the loci with at least n_proj called allele copies in every population are used. A folded spectrum is folded on the total number of copies across all populations (as in dadi).

Usage

nwise_pop_sfs(.x, ...)

## S3 method for class 'gen_tbl'
nwise_pop_sfs(
  .x,
  n_proj = NULL,
  folded = TRUE,
  ancestral = NULL,
  type = c("tidy", "array"),
  fastest = c("last", "first"),
  max_cells = 1e+07,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

Arguments

Details

Populations are in the order of the groups (see dplyr::group_keys()): to change their order (which matters for dadi and fastsimcoal2), change the order of the levels of the grouping variable (e.g. with factor()); to use a subset of the populations, filter() the gen_tibble (empty groups are dropped). The spectrum has prod(n_proj + 1) cells, which grows quickly with the number of populations.

Value

If type = "tidy", a tibble with one row per cell of the joint spectrum, with one column per population (named with the group labels) giving its number of copies of the derived (or minor) allele (see allele_count in pop_sfs()), n_sites (the projected number of sites) and masked (whether the cell is masked, see pop_sfs()). The projection of each population is stored in the attribute n_proj. If type = "array", an sfs object (see pop_sfs()), an array with one dimension per population.

See Also

pop_sfs(), pairwise_pop_sfs(), pop_sfs_preview()

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

example_gt %>% nwise_pop_sfs(n_proj = 2)

# the first population changing fastest
example_gt %>% nwise_pop_sfs(n_proj = 2, fastest = "first")

Compute the Pairwise Allele Sharing Matrix for a gen_tibble object

Description

This function computes the Allele Sharing matrix. Estimates Allele Sharing (equivalent to the quantity estimated by hierfstat::matching()) between pairs of individuals (for each locus, gives 1 if the two individuals are homozygous for the same allele, 0 if they are homozygous for a different allele, and 1/2 if at least one individual is heterozygous. Matching is the average of these 0, 1/2 and 1s)

Usage

pairwise_allele_sharing(
  x,
  as_matrix = FALSE,
  block_size = bigstatsr::block_size(nrow(x)),
  n_cores = bigstatsr::nb_cores()
)

Arguments

Details

Imputed genotypes are used if x is set to use them (see gt_set_imputed()), and treated as missing otherwise.

Value

a matrix of allele sharing between all pairs of individuals

See Also

hierfstat::matching()

Examples

example_gt <- load_example_gt("gen_tbl")

# Compute allele sharing between individuals
example_gt %>% pairwise_allele_sharing(as_matrix = FALSE)

# Alternatively, return as a tibble
example_gt %>% pairwise_allele_sharing(as_matrix = TRUE)

Compute the Genomic Relationship Matrix for a gen_tibble object

Description

This function computes the Genomic Relationship Matrix (GRM). This is estimated by computing the pairwise kinship coefficients (coancestries) between all pairs of individuals from a matrix of Allele Sharing following the approach of Weir and Goudet 2017 based on beta estimators).

Usage

pairwise_grm(
  x,
  allele_sharing_mat = NULL,
  block_size = bigstatsr::block_size(nrow(x)),
  n_cores = bigstatsr::nb_cores()
)

Arguments

Details

The GRM is twice the coancestry matrix (e.g. as estimated by hierfstat::beta.dosage() with inb=FALSE).

Value

a matrix of GR between all pairs of individuals

See Also

hierfstat::beta.dosage()

Examples

example_gt <- load_example_gt("gen_tbl")

# Compute the GRM from the allele sharing matrix
example_gt %>% pairwise_grm()

# To calculate using a precomputed allele sharing matrix, use:
allele_sharing <- example_gt %>% pairwise_allele_sharing(as_matrix = TRUE)
example_gt %>% pairwise_grm(allele_sharing_mat = allele_sharing)

Compute the Identity by State Matrix for a gen_tibble object

Description

This function computes the IBS matrix.

Usage

pairwise_ibs(
  x,
  as_matrix = FALSE,
  type = c("proportion", "adjusted_counts", "raw_counts"),
  block_size = bigstatsr::block_size(nrow(x)),
  n_cores = bigstatsr::nb_cores()
)

Arguments

Details

Note that monomorphic sites are currently considered (as in PLINK). Remove monomorphic sites before running pairwise_ibs() if this is a concern.

Imputed genotypes are used if x is set to use them (see gt_set_imputed()), and treated as missing otherwise.

Missing genotypes are handled as in PLINK with the flat-missing modifier (⁠--distance ibs flat-missing⁠ for "proportion", ⁠--distance flat-missing⁠ for "adjusted_counts"): each pair of individuals is compared only on the loci genotyped in both, and all loci have the same weight. By default (without flat-missing), PLINK 1.9 instead divides each distance by one minus the sum of the average contributions to the distance of the loci missing in the pair, which depend on their allele frequencies. This default correction is not available in pairwise_ibs(), so, when genotypes are missing, its values are only comparable to those of PLINK run with flat-missing. Without missing genotypes, both corrections give the same values.

Value

for "proportion" and "adjusted_counts", a tidy tibble with one row per pair of individuals (columns item1, item2 and value), or a square matrix if as_matrix = TRUE; for "raw_counts", a list of two bigstatsr::FBM matrices, ibs (the counts of alleles identical by state) and valid_n (the counts of non-missing alleles, i.e. twice the number of loci genotyped in both individuals).

Examples

example_gt <- load_example_gt("gen_tbl")

pairwise_ibs(example_gt, type = "proportion")

# Alternatively, return a matrix
pairwise_ibs(example_gt, type = "proportion", as_matrix = TRUE)

# Adjust block_size
pairwise_ibs(example_gt, block_size = 2)

# Change type
pairwise_ibs(example_gt, type = "adjusted_counts")
pairwise_ibs(example_gt, type = "raw_counts")

Compute the KING-robust Matrix for a gen_tibble object

Description

This function computes the KING-robust estimator of kinship, reimplementing the KING algorithm of Manichaikul et al. (2010).

Usage

pairwise_king(
  x,
  as_matrix = FALSE,
  block_size = bigstatsr::block_size(nrow(x)),
  n_cores = bigstatsr::nb_cores()
)

Arguments

Details

Imputed genotypes are used if x is set to use them (see gt_set_imputed()), and treated as missing otherwise.

Value

a square symmetrical matrix of relationship coefficients between individuals if as_matrix is TRUE, or a tidied tibble of coefficients if as_matrix is FALSE.

References

Manichaikul, A. et al. (2010) Robust relationship inference in genome-wide association studies. Bioinformatics, 26(22), 2867–2873. https://doi.org/10.1093/bioinformatics/btq559.

Note that monomorphic sites are currently considered. Remove monomorphic sites before running pairwise_king if this is a concern.

Examples

example_gt <- load_example_gt("gen_tbl")

# Compute the KING-robust matrix
pairwise_king(example_gt, as_matrix = TRUE)

# Or return a tidy tibble
pairwise_king(example_gt, as_matrix = FALSE)

# Adjust block_size
pairwise_king(example_gt, block_size = 2)

Compute pairwise population Fst

Description

This function computes pairwise Fst. The following methods are implemented:

• 'Hudson': Hudson's formulation, as derived in Bhatia et al (2013) for diploids. This is the only method that can also be used with pseudohaploid data.

• 'Nei87' : Fst according to Nei (1987) - includes the correction for heterozygosity when computing Ht (it uses the same formulation as in hierfstat::pairwise.neifst()),

• 'WC84' : Weir and Cockerham (1984), correcting for missing data (it uses the same formulation as in hierfstat::pairwise.WCfst()).

Usage

pairwise_pop_fst(.x, ...)

## S3 method for class 'gen_tbl'
pairwise_pop_fst(
  .x,
  type = c("tidy", "pairwise"),
  by_locus = FALSE,
  by_locus_type = c("tidy", "matrix", "list"),
  method = c("Hudson", "Nei87", "WC84"),
  return_num_dem = FALSE,
  n_cores = bigstatsr::nb_cores(),
  ...
)

Arguments

Details

For all formulae, the genome wide estimate is obtained by taking the ratio of the mean numerators and denominators over all relevant SNPs.

Loci with no genotypes in one of the two populations of a pair (e.g. loci that failed in all the individuals of a population) carry no information on the differentiation between them, so they are dropped for that pair with all methods: their Fst by locus is NaN, and they do not contribute to the mean numerator and denominator of the genome wide estimate.

Difference from hierfstat for 'WC84': hierfstat::pairwise.WCfst() (and hierfstat::wc()) keep these loci. Their among-population variance component (\sigma^2_a, the numerator) is set to 0, while their within-population components (\sigma^2_b and \sigma^2_w, estimated from the one population with data) still add to the denominator. So hierfstat gives these loci an Fst of 0 rather than NaN, and its genome wide estimate is pulled towards 0 (i.e. it is usually lower than the one computed here). The two estimates are identical when every locus has data in both populations, or if such loci are removed before using hierfstat. This applies to every function that computes Fst with 'WC84' (e.g. windows_pairwise_pop_fst(), jackknife_pairwise_pop_fst() and nwise_pop_pbs()). 'Nei87' gives the same results as hierfstat::pairwise.neifst(), which also drops these loci.

Value

if type=tidy, a tibble of genome-wide pairwise Fst values with each pairwise combination as a row if "by_locus=FALSE", else a list including the tibble of genome-wide values as well as a matrix with pairwise Fst by locus with loci as rows and and pairwise combinations as columns. If type=pairwise, a matrix of genome-wide pairwise Fst values is returned.

References

Bhatia G, Patterson N, Sankararaman S, Price AL. (2013) Estimating and Interpreting FST: The Impact of Rare Variants. Genome Research, 23(9):1514–1521.

Nei, M. (1987) Molecular Evolutionary Genetics. Columbia University Press

Weir, B. S., & Cockerham, C. C. (1984). Estimating F-statistics for the analysis of population structure. Evolution, 38(6): 1358–1370.

See Also

hierfstat::pairwise.neifst()

Examples

example_gt <- load_example_gt("gen_tbl")

# For a basic global pairwise Fst calculation:
example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Nei87")

# With a pairwise matrix:
example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Nei87", type = "pairwise")

# To calculate Fst by locus:
example_gt %>%
  group_by(population) %>%
  pairwise_pop_fst(method = "Hudson", by_locus = TRUE)

Compute the joint site frequency spectrum of each pair of populations

Description

Compute the two-dimensional (joint) site frequency spectrum (SFS) of each pair of populations, i.e. the number of loci with i copies of the derived (or minor) allele in the first population and j in the second. Missing genotypes are dealt with as in easySFS, by projecting each population down to n_proj allele copies (see pop_sfs() for details). Each spectrum only uses the loci that have at least n_proj called allele copies in both its populations, so its marginal spectra are not, in general, equal to the spectra computed by pop_sfs(). A folded spectrum is folded on the total number of copies across the two populations (as in dadi), so that it counts copies of the allele that is the minor one across the pair.

Usage

pairwise_pop_sfs(.x, ...)

## S3 method for class 'gen_tbl'
pairwise_pop_sfs(
  .x,
  n_proj = NULL,
  folded = TRUE,
  ancestral = NULL,
  type = c("tidy", "array"),
  fastest = c("last", "first"),
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

Arguments

Value

If type = "tidy", a tibble with one row per cell of the joint spectrum of each pair of populations, with the columns:

• the groups of the pair (e.g. population_1 and population_2; one column per grouping variable for each population of the pair, if the data are grouped by more than one variable);

• allele_count_1 and allele_count_2: the number of copies of the derived (or minor) allele in the first and second population (see allele_count in pop_sfs());

• n_sites: the (projected) number of sites;

• masked: whether the cell is masked (see pop_sfs()).

The projection of each population is stored in the attribute n_proj. The tibble has class pairwise_sfs_tbl, and can be plotted with autoplot(). If type = "array", a list of sfs objects (one per pair, named with the labels of the two groups pasted with "."), each an array with the first population in the rows and the second in the columns (see pop_sfs()). Pairs are in the order given by utils::combn().

See Also

pop_sfs(), nwise_pop_sfs(), pop_sfs_preview()

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

example_gt %>% pairwise_pop_sfs(n_proj = 2)

# as arrays
example_gt %>% pairwise_pop_sfs(n_proj = 2, type = "array")

Compute population specific FIS

Description

This function computes population specific FIS, using either the approach of Nei 1987 (with an algorithm equivalent to the one used by hierfstat::basic.stats()) or of Weir and Goudet 2017 (with an algorithm equivalent to the one used by hierfstat::fis.dosage()).

Usage

pop_fis(
  .x,
  method = c("Nei87", "WG17"),
  by_locus = FALSE,
  include_global = FALSE,
  allele_sharing_mat = NULL,
  n_cores = bigstatsr::nb_cores()
)

Arguments

Value

a vector of population specific fis (plus the global value if include_global=TRUE)

References

Nei M. (1987) Molecular Evolutionary Genetics. Columbia University Press.

Weir, BS and Goudet J (2017) A Unified Characterization of Population Structure and Relatedness. Genetics (2017) 206:2085

See Also

hierfstat::basic.stats() hierfstat::fis.dosage()

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# Compute FIS using Nei87:
example_gt %>% pop_fis(method = "Nei87")

# Compute FIS using WG17:
example_gt %>% pop_fis(method = "WG17")

# To include the global FIS, set include_global = TRUE:
example_gt %>% pop_fis(method = "Nei87", include_global = TRUE)

# To return FIS by locus, set by_locus = TRUE:
example_gt %>% pop_fis(method = "Nei87", by_locus = TRUE)

# To calculate from a pre-computed allele sharing matrix:
allele_sharing_mat <- pairwise_allele_sharing(example_gt, as_matrix = TRUE)
example_gt %>% pop_fis(
  method = "WG17", allele_sharing_mat =
    allele_sharing_mat
)

Compute population specific Fst

Description

This function computes population specific Fst, using the approach in Weir and Goudet 2017 (as computed by hierfstat::fst.dosage()).

Usage

pop_fst(.x, include_global = FALSE, allele_sharing_mat = NULL)

Arguments

Value

a vector of population specific Fst (plus the global value if include_global=TRUE)

References

Weir, BS and Goudet J (2017) A Unified Characterization of Population Structure and Relatedness. Genetics (2017) 206:2085

See Also

hierfstat::fst.dosage()

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# Compute FIS using Nei87
example_gt %>% pop_fst()

# To include the global Fst, set include_global = TRUE
example_gt %>% pop_fst(include_global = TRUE)

# To calculate from a pre-computed allele sharing matrix:
allele_sharing_mat <- pairwise_allele_sharing(example_gt, as_matrix = TRUE)
example_gt %>% pop_fst(allele_sharing_mat = allele_sharing_mat)

Compute basic population global statistics

Description

This function computes basic population global statistics, following the notation in Nei 1987 (which in turn is based on Nei and Chesser 1983):

• observed heterozygosity ( \hat{h}_o, column header Ho)

• expected heterozygosity, also known as gene diversity ( \hat{h}_s, Hs)

• total heterozygosity ( \hat{h}_t, Ht)

• genetic differentiation between subpopulations (D_{st}, Dst)

• corrected total population diversity (h'_t, Htp)

• corrected genetic differentiation between subpopulations (D'_{st}, Dstp)

• \hat{F}_{ST} (column header, Fst)

• corrected \hat{F'}_{ST} (column header Fstp)

• \hat{F}_{IS} (column header, Fis)

• Jost's \hat{D} (column header, Dest)

Usage

pop_global_stats(.x, by_locus = FALSE, n_cores = bigstatsr::nb_cores())

Arguments

Details

We use the notation of Nei 1987. That notation was for loci with m alleles, but in our case we only have two alleles, so m=2.

• Within population observed heterozygosity \hat{h}_o for a locus with m alleles is defined as:
  \hat{h}_o= 1-\sum_{k=1}^{s} \sum_{i=1}^{m} \hat{X}_{kii}/s
  where
  \hat{X}_{kii} represents the proportion of homozygote i in the sample for the kth population and
  s the number of populations,
  following equation 7.38 in Nei(1987) on pp.164.
  

• Within population expected heterozygosity (gene diversity) \hat{h}_s for a locus with m alleles is defined as:
  \hat{h}_s=(\tilde{n}/(\tilde{n}-1))[1-\sum_{i=1}^{m}\bar{\hat{x}_i^2}-\hat{h}_o/2\tilde{n}]
  #nolint where
  \tilde{n}=s/\sum_k 1/n_k (i.e the harmonic mean of n_k) and
  \bar{\hat{x}_i^2}=\sum_k \hat{x}_{ki}^2/s
  following equation 7.39 in Nei(1987) on pp.164.

• Total heterozygosity (total gene diversity) \hat{h}_t for a locus with m alleles is defined as:
  \hat{h}_t = 1-\sum_{i=1}^{m} \bar{\hat{x}_i^2} + \hat{h}_s/(\tilde{n}s) - \hat{h}_o/(2\tilde{n}s)
  where
  \hat{x}_i=\sum_k \hat{x}_{ki}/s
  following equation 7.40 in Nei(1987) on pp.164.
  

• The amount of gene diversity among samples D_{ST} is defined as:
  D_{ST} = \hat{h}_t - \hat{h}_s
  following the equation provided in the text at the top of page 165 in Nei(1987).

• The corrected amount of gene diversity among samples D'_{ST} is defined as:
  D'_{ST} = (s/(s-1))D'_{ST}
  following the equation provided in the text at the top of page 165 in Nei(1987).

• Total corrected heterozygosity (total gene diversity) \hat{h}_t is defined as:
  \hat{h'}_t = \hat{h}_s + D'_{ST}
  following the equation provided in the text at the top of page 165 in Nei(1987).

• \hat{F}_{IS} is defined as:
  \hat{F}_{IS} = 1 - \hat{h}_o/\hat{h}_s
  following equation 7.41 in Nei(1987) on pp.164.
  

• \hat{F}_{ST} is defined as:
  \hat{F}_{ST} = 1 - \hat{h}_s/\hat{h}_t = D_{ST}/\hat{h}_t
  following equation 7.43 in Nei(1987) on pp.165.
  

• \hat{F'}_{ST} is defined as:
  \hat{F'}_{ST} = D'_{ST}/\hat{h'}_t
  following the explanation provided in the text at the top of page 165 in Nei(1987).

• Jost's \hat{D} is defined as:
  \hat{D} = (s/(s-1))((\hat{h'}_t-\hat{h}_s)/(1-\hat{h}_s))
  as defined by Jost(2008)

  All these statistics are first computed by locus, and then averaged across loci (including any monomorphic locus) to obtain genome-wide values. The function uses the same algorithm as hierfstat::basic.stats() but is optimized for speed and memory usage.

Value

a tibble of population statistics, with populations as rows and statistics as columns

References

Nei M, Chesser R (1983) Estimation of fixation indexes and gene diversities. Annals of Human Genetics, 47, 253-259.

Nei M. (1987) Molecular Evolutionary Genetics. Columbia University Press, pp. 164-165.

Jost L (2008) GST and its relatives do not measure differentiation. Molecular Ecology, 17, 4015-4026.

See Also

hierfstat::basic.stats()

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# Compute population global statistics
example_gt %>% pop_global_stats()

# To return by locus, set by_locus = TRUE
example_gt %>% pop_global_stats(by_locus = TRUE)

Compute the population expected heterozygosity

Description

This function computes expected population heterozygosity (also referred to as gene diversity, to avoid the potentially misleading use of the term "expected" in this context), using the formula of Nei (1987).

Usage

pop_het_exp(
  .x,
  by_locus = FALSE,
  include_global = FALSE,
  n_cores = bigstatsr::nb_cores()
)

pop_gene_div(
  .x,
  by_locus = FALSE,
  include_global = FALSE,
  n_cores = bigstatsr::nb_cores()
)

Arguments

Details

Within population expected heterozygosity (gene diversity) \hat{h}_s for a locus with m alleles is defined as:
\hat{h}_s=\tilde{n}/(\tilde{n}-1)[1-\sum_{i}^{m}\bar{\hat{x}_i^2}-\hat{h}_o/2\tilde{n}]
#nolint

where
\tilde{n}=s/\sum_k 1/n_k (i.e the harmonic mean of n_k) and
\bar{\hat{x}_i^2}=\sum_k \hat{x}_{ki}^2/s
following equation 7.39 in Nei(1987) on pp.164. In our specific case, there are only two alleles, so m=2. \hat{h}_s at the genome level for each population is simply the mean of the locus estimates for each population.

Value

a vector of mean population observed heterozygosities (if by_locus=FALSE), or a matrix of estimates by locus (rows are loci, columns are populations, by_locus=TRUE)

References

Nei M. (1987) Molecular Evolutionary Genetics. Columbia University Press

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# Compute expected heterozygosity
example_gt %>% pop_het_exp()

# To include the global expected heterozygosity, set include_global = TRUE
example_gt %>% pop_het_exp(include_global = TRUE)

# To return by locus, set by_locus = TRUE
example_gt %>% pop_het_exp(by_locus = TRUE)

Compute the population observed heterozygosity

Description

This function computes population heterozygosity, using the formula of Nei (1987).

Usage

pop_het_obs(
  .x,
  by_locus = FALSE,
  include_global = FALSE,
  n_cores = bigstatsr::nb_cores()
)

Arguments

Details

Within population observed heterozygosity \hat{h}_o for a locus with m alleles is defined as:
\hat{h}_o= 1-\sum_{k=1}^{s} \sum_{i=1}^{m} \hat{X}_{kii}/s
where
\hat{X}_{kii} represents the proportion of homozygote i in the sample for the kth population and
s the number of populations,
following equation 7.38 in Nei(1987) on pp.164. In our specific case, there are only two alleles, so m=2. For population specific estimates, the sum is done over a single value of k. \hat{h}_o at the genome level is simply the mean of the locus estimates.

This statistic is ploidy-agnostic: an individual is heterozygous at a locus whenever its dosage is neither zero nor its own ploidy, so pop_het_obs() works for diploid, polyploid, and populations that mix individuals of different ploidy. The genome-wide estimate when include_global = TRUE is the unweighted mean of the population estimates for polyploid/mixed-ploidy data; for diploid and pseudohaploid data it is computed via pop_global_stats(), for consistency with previous releases.

Value

a vector of mean population observed heterozygosities (if by_locus=FALSE), or a matrix of estimates by locus (rows are loci, columns are populations, by_locus=TRUE)

References

Nei M. (1987) Molecular Evolutionary Genetics. Columbia University Press

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# Compute expected heterozygosity
example_gt %>% pop_het_obs()

# To include the global expected heterozygosity, set include_global = TRUE
example_gt %>% pop_het_obs(include_global = TRUE)

# To return by locus, set by_locus = TRUE
example_gt %>% pop_het_obs(by_locus = TRUE)

Compute the site frequency spectrum of each population

Description

Compute the one-dimensional site frequency spectrum (SFS) of each population, i.e. the number of loci with 0, 1, 2, ... copies of the derived (or minor) allele. Missing genotypes are dealt with as in easySFS (which relies on dadi), so that the spectra can be used for demographic inference with dadi, moments or fastsimcoal2. See pairwise_pop_sfs() for the joint spectra of pairs of populations, and nwise_pop_sfs() for the joint spectrum of all populations.

Usage

pop_sfs(.x, ...)

## S3 method for class 'tbl_df'
pop_sfs(
  .x,
  n_proj = NULL,
  folded = TRUE,
  ancestral = NULL,
  type = c("tidy", "array"),
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

## S3 method for class 'vctrs_bigSNP'
pop_sfs(
  .x,
  n_proj = NULL,
  folded = TRUE,
  ancestral = NULL,
  type = c("tidy", "array"),
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(length(.x), 1),
  ...
)

## S3 method for class 'grouped_gen_tbl'
pop_sfs(
  .x,
  n_proj = NULL,
  folded = TRUE,
  ancestral = NULL,
  type = c("tidy", "array"),
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

Arguments

Details

Missing data. Because of missing genotypes, the number of allele copies called at a locus varies from locus to locus. Each population is therefore projected down to a fixed number of allele copies, n_proj (Marth et al. 2004; Gutenkunst et al. 2009): a locus with n called allele copies, k of which are derived, contributes to cell j of the spectrum the hypergeometric probability of drawing j derived copies when sampling n_proj of its n copies without replacement. Loci with fewer than n_proj called copies are dropped. Each locus is thus spread over several cells, and the spectrum contains fractional numbers of sites. This is the approach of easySFS (with the -a option, i.e. using all loci).

Choosing the projection. Projecting to more allele copies gives more resolution (more segregating sites, especially rare variants), but drops more loci with missing data. Following the advice of the dadi manual and easySFS, use pop_sfs_preview() to compute the number of segregating sites for each possible projection, and choose a value that keeps most segregating sites. If n_proj is not given, each population is projected to the number of allele copies that maximises its number of segregating sites (the best projection in pop_sfs_preview()).

Folding and polarisation. Unless folded = FALSE, the spectrum is folded (i.e. it counts copies of the minor allele), as the ancestral allele is not usually known. An unfolded spectrum counts copies of the derived allele. By default (ancestral = NULL), the reference allele is taken to be ancestral, as in easySFS, so that the spectrum counts copies of the alternate allele: this is only meaningful if the alleles of the gen_tibble are polarised (note that some operations, such as merging datasets with rbind(), can swap the reference and alternate alleles). Otherwise, ancestral gives the ancestral allele of each locus (e.g. from an outgroup, or the AA field of a vcf), either as a column of the loci table or as a vector: loci whose ancestral allele is the alternate allele are polarised (copies of the reference allele are counted), and loci whose ancestral allele is missing (NA) or is neither of their two alleles are dropped, as in dadi. Alleles are compared ignoring case. Polarisation does not change a folded spectrum, so ancestral is ignored (with a warning) when folded = TRUE.

Masked cells. As in dadi, the corners of the spectrum (sites monomorphic in the projected sample) are masked, and so are the cells that were folded out in a folded spectrum. They are reported in the masked column (or the mask attribute), but keep their values.

Each spectrum only uses the loci that have enough called allele copies in its own populations. Diploid, polyploid, mixed-ploidy and pseudohaploid data are supported: each individual contributes as many allele copies as its ploidy (one for pseudohaploids). The spectrum is computed from the genotypes currently in use: if the gen_tibble uses imputed genotypes (see gt_set_imputed()), these are treated as called genotypes, so you will usually want to use gt_set_imputed(x, FALSE) first.

Note that easySFS samples a single SNP per CHROM (i.e. per RAD locus) by default; to do the same, thin the loci (e.g. with select_loci()) before computing the spectra.

Value

If type = "tidy", a tibble with one row per cell of the spectrum of each population, with the columns:

• group (or one column per grouping variable, if the data are grouped by more than one variable), absent for ungrouped data;

• allele_count: the number of copies of the derived allele (for an unfolded spectrum: the alternate allele, unless the loci are polarised with ancestral) or of the minor allele (for a folded spectrum);

• n_sites: the (projected) number of sites;

• masked: whether the cell is masked (see details).

The projection of each population is stored in the attribute n_proj. The tibble has class pop_sfs_tbl, and can be plotted with autoplot(). If type = "array", a list of sfs objects (one per population, named by the group labels), or a single sfs object for ungrouped data. An sfs object is an array of the numbers of sites, with attributes giving its mask, the population labels (pop_ids), n_proj, whether it is folded, and the number of loci used (n_loci), dropped for lack of called allele copies (n_loci_dropped), and dropped because their ancestral allele is unknown (n_loci_unpolarised, only with ancestral). Use tidy() to turn it into a tibble.

References

Marth GT, Czabarka E, Murvai J, Sherry ST (2004) The allele frequency spectrum in genome-wide human variation data reveals signals of differential demographic history in three large world populations. Genetics 166: 351-372.

Gutenkunst RN, Hernandez RD, Williamson SH, Bustamante CD (2009) Inferring the joint demographic history of multiple populations from multidimensional SNP frequency data. PLoS Genetics 5: e1000695.

Overcast I. easySFS. https://github.com/isaacovercast/easySFS

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# Preview the number of segregating sites for each projection
example_gt %>% pop_sfs_preview()

# Compute the folded spectra, projecting each population to 2 copies
example_gt %>% pop_sfs(n_proj = 2)

# Project each population to a different number of copies
example_gt %>% pop_sfs(n_proj = c(pop1 = 4, pop2 = 2, pop3 = 2))

# Unfolded spectra, as arrays
example_gt %>% pop_sfs(n_proj = 2, folded = FALSE, type = "array")

# Unfolded spectra, polarised with the ancestral allele of each locus (the
# locus with an unknown ancestral allele is dropped)
show_loci(example_gt)
example_gt %>%
  pop_sfs(
    n_proj = 2, folded = FALSE,
    ancestral = c("T", "C", "C", "G", NA, "T")
  )

Preview the projections of the site frequency spectrum

Description

For each population, compute the number of segregating sites in its site frequency spectrum projected to each possible number of allele copies (from 2 to the number of allele copies of the population), as the --preview option of easySFS. This is used to choose n_proj in pop_sfs(), pairwise_pop_sfs() and nwise_pop_sfs(): projecting to more allele copies gives more resolution (more segregating sites, especially rare variants), but drops more loci with missing data. The dadi manual recommends maximising the number of segregating sites, while keeping a reasonable sample size.

Usage

pop_sfs_preview(.x, ...)

## S3 method for class 'tbl_df'
pop_sfs_preview(
  .x,
  ancestral = NULL,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

## S3 method for class 'vctrs_bigSNP'
pop_sfs_preview(
  .x,
  ancestral = NULL,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(length(.x), 1),
  ...
)

## S3 method for class 'grouped_gen_tbl'
pop_sfs_preview(
  .x,
  ancestral = NULL,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

Arguments

Details

The number of segregating sites for a projection to p copies is the expected number of loci that are polymorphic in a sample of p copies drawn without replacement from the called copies of each locus (loci with fewer than p called copies are dropped), i.e. the sum of the cells of the projected spectrum that are not corners. It is the same for the folded and unfolded spectrum. If unfolded spectra are polarised with ancestral (see pop_sfs()), give the same ancestral here, so that the loci without a known ancestral allele are left out of the preview, as they are of the spectra.

Value

a tibble with columns group (or one column per grouping variable, if the data are grouped by more than one variable; absent for ungrouped data), n_proj (the number of allele copies of the projection), n_seg_sites (the number of segregating sites, not rounded) and n_loci (the number of loci with at least n_proj called allele copies), and best, which is TRUE for the projection that maximises the number of segregating sites (the largest one in case of ties), used by default by pop_sfs(), pairwise_pop_sfs() and nwise_pop_sfs(). The tibble has class sfs_preview_tbl, and can be plotted with autoplot().

See Also

pop_sfs()

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

example_gt %>% pop_sfs_preview()

# the projection chosen for each population
example_gt %>%
  pop_sfs_preview() %>%
  dplyr::filter(best)

# plot the number of segregating sites for each projection
example_gt %>%
  pop_sfs_preview() %>%
  autoplot()

Estimate Tajima's D for the whole genome

Description

Note that Tajima's D estimates from data that have been filtered or ascertained can be difficult to interpret. This function should ideally be used on sequence data prior to filtering.

Usage

pop_tajimas_d(.x, ...)

## S3 method for class 'tbl_df'
pop_tajimas_d(
  .x,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

## S3 method for class 'vctrs_bigSNP'
pop_tajimas_d(
  .x,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(length(.x), 1),
  ...
)

## S3 method for class 'grouped_gen_tbl'
pop_tajimas_d(
  .x,
  n_cores = bigstatsr::nb_cores(),
  block_size = bigstatsr::block_size(nrow(.x), 1),
  ...
)

Arguments

Details

Nucleotide diversity is estimated at each locus with loci_pi(), so this function works with diploid, polyploid, mixed-ploidy and pseudohaploid data.

Missing genotypes are accounted for following Bailey et al. (2025), as implemented in pixy: the sample size n_x at each locus x is the number of non-missing allele copies at that locus (i.e. the sum of the ploidy of the genotyped individuals, with pseudohaploids contributing one copy each). The two estimators of \theta are then:

• \hat{\theta}_\pi = \sum_x \pi_x, where \pi_x is the nucleotide diversity at locus x, estimated from its n_x allele copies (loci with fewer than 2 allele copies contribute 0);

• \hat{\theta}_W = \sum_x 1 / a_{n_x}, summed over the segregating sites, where a_{n_x} = \sum_{i=1}^{n_x - 1} 1/i (i.e. each segregating site is weighted by its own sample size).

The standard deviation of \hat{\theta}_\pi - \hat{\theta}_W is computed with the formula of Tajima (1989), using the number of segregating sites S and the mean sample size across loci with at least one non-missing allele copy (rounded to the nearest integer, as the formula is only defined for integer sample sizes). When there are no missing data, all loci have the same sample size and the statistic is identical to the classic Tajima's D. Loci with no genotyped individuals are ignored. Tajima's D is NA when there are no segregating sites, or when the mean sample size is smaller than 2.

Value

A single numeric value (Tajima's D D) for the whole data set, NA when the statistic is not defined. For grouped data a list of Tajima's D D values (one per group) is returned.

References

Tajima F (1989) Statistical method for testing the neutral mutation hypothesis by DNA polymorphism. Genetics 123:585-595.

Bailey N, Stevison L, Samuk K (2025) Correcting for bias in estimates of \theta_W and Tajima's D from missing data in next-generation sequencing. Molecular Ecology Resources, e14104. doi:10.1111/1755-0998.14104

Examples

example_gt <- load_example_gt("grouped_gen_tbl")

# Compute Tajima's D
example_gt %>% pop_tajimas_d()

Predict scores of a PCA

Description

Predict the PCA scores for a gt_pca, either for the original data or projecting new data.

Usage

## S3 method for class 'gt_pca'
predict(
  object,
  new_data = NULL,
  project_method = c("none", "simple", "OADP", "least_squares"),
  lsq_pcs = c(1, 2),
  block_size = NULL,
  n_cores = 1,
  as_matrix = TRUE,
  ...
)

Arguments

Value

a matrix of predictions (in line with predict using a prcomp object) or a tibble, with samples as rows and components as columns. The number of components depends on how many were estimated in the gt_pca object.

References

Zhang et al (2020). Fast and robust ancestry prediction using principal component analysis 36(11): 3439–3446.

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Remove monomorphic loci and impute
lobsters <- lobsters %>% select_loci_if(loci_maf(genotypes) > 0)
lobsters <- gt_impute_simple(lobsters, method = "mode")

# Subset into two datasets: one original and one to predict
original_lobsters <- lobsters[c(1:150), ]
new_lobsters <- lobsters[c(151:176), ]

# Create PCA object
pca <- gt_pca_partialSVD(original_lobsters)

# Predict
predict(pca, new_data = new_lobsters, project_method = "simple")

# Predict with OADP
predict(pca, new_data = new_lobsters, project_method = "OADP")

# Predict with least squares
predict(pca,
  new_data = new_lobsters,
  project_method = "least_squares", lsq_pcs = c(1, 2, 3)
)

# Return a tibble
predict(pca, new_data = new_lobsters, as_matrix = FALSE)

# Adjust block.size
predict(pca, new_data = new_lobsters, block_size = 10)

Print method for matrix_int_names

Description

Print method for matrix_int_names

Usage

## S3 method for class 'matrix_int_names'
print(x, ..., n = 10L)

Arguments

Value

Invisibly returns the original object

See Also

Other matrix_int_names_functions: [.matrix_int_names(), as.data.frame.matrix_int_names(), col_names(), matrix_int_names(), row_names()

Examples

my_mat <- matrix_int_names(matrix(1:6, nrow = 2),
  row_names = c(10L, 20L),
  col_names = c(100L, 200L, 300L)
)
print(my_mat)
print(my_mat, n = 2)

Convert a standard matrix to a q_matrix object

Description

Takes a single Q matrix that exists as either a matrix or a data frame and returns a q_matrix object.

Usage

q_matrix(x)

Arguments

Value

A q_matrix object

Examples

# Read in a single .Q file
q_mat <- read.table(system.file("extdata", "anolis", "anolis_ld_run1.3.Q",
  package = "tidypopgen"
))
class(q_mat)

# Convert to a Q matrix object
q_mat <- q_matrix(q_mat)
class(q_mat)

Create a Quality Control report for individuals

Description

Return QC information to assess individuals (Observed heterozygosity and missingness).

Usage

qc_report_indiv(.x, ...)

## S3 method for class 'gen_tbl'
qc_report_indiv(.x, kings_threshold = NULL, ...)

## S3 method for class 'grouped_gen_tbl'
qc_report_indiv(.x, kings_threshold = NULL, ...)

Arguments

Details

Providing the parameter kings_threshold will return two additional columns, 'id' containing the ID of individuals, and 'to_keep' a logical vector describing whether the individual should be removed to retain the largest possible set of individuals with no relationships above the threshold. The calculated pairwise KING relationship matrix is also returned as an attribute of 'to_keep'. The kings_threshold parameter can be either a numeric KING kinship coefficient or a string of either "first" or "second", to remove any first degree or second degree relationships from the dataset. This second option is similar to using –unrelated –degree 1 or –unrelated –degree 2 in KING. For pseudohaploid data, only missingness and ploidy are reported.

Value

If no kings_threshold is provided, a tibble with 2 elements: het_obs and missingness. If kings_threshold is provided, a tibble with 4 elements: het_obs, missingness, id and to_keep. For pseudohaploid data, a tibble with ploidy and missingness.

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
example_gt <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Get QC report for individuals
example_gt %>% qc_report_indiv()

# Get QC report with kinship filtering
example_gt %>% qc_report_indiv(kings_threshold = "first")

Create a Quality Control report for loci

Description

Return QC information to assess loci (MAF, missingness and HWE test). For pseudohaploid data, HWE test is not calculated.

Usage

qc_report_loci(.x, ...)

## S3 method for class 'gen_tbl'
qc_report_loci(.x, ...)

## S3 method for class 'grouped_gen_tbl'
qc_report_loci(.x, ...)

Arguments

Value

either a tibble with 3 elements (maf, missingness and hwe_p). For pseudohaploid data, a tibble with 2 elements (maf and missingness).

Examples

# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
example_gt <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Get a QC report for the loci
example_gt %>% qc_report_loci()

# Group by population to calculate HWE within populations
example_gt <- example_gt %>% group_by(population)
example_gt %>% qc_report_loci()

Combine two gen_tibbles

Description

This function combined two gen_tibbles. It subsets the loci and swaps ref and alt alleles to make the two datasets compatible. The first object is used as a "reference" , and SNPs in the other dataset will be flipped and/or alleles swapped as needed. SNPs that have different alleles in the two datasets (i.e. triallelic) will also be dropped. There are also options (NOT default) to attempt strand flipping to match alleles (often needed in human datasets from different SNP chips), and remove ambiguous alleles (C/G and A/T) where the correct strand can not be guessed.

Usage

## S3 method for class 'gen_tbl'
rbind(
  ...,
  as_is = FALSE,
  flip_strand = FALSE,
  use_position = FALSE,
  quiet = FALSE,
  backingfile = NULL
)

Arguments

Details

Loci are matched by name, unless use_position = TRUE. If the loci names of either gen_tibble are integers (as created by gen_tibble(names_as_int = TRUE)), loci are always matched by chromosome and position (see use_position).

rbind differs from merging data with plink, which swaps the order of allele1 and allele2 according to minor allele frequency when merging datasets. rbind flips and/or swaps alleles according to the reference dataset, not according to allele frequency.

Value

a gen_tibble with the merged data.

Examples

example_gt <- load_example_gt("gen_tbl")

# Create a second gen_tibble to merge
test_indiv_meta <- data.frame(
  id = c("x", "y", "z"),
  population = c("pop1", "pop1", "pop2")
)
test_genotypes <- rbind(
  c(1, 1, 0, 1, 1, 0),
  c(2, 1, 0, 0, 0, 0),
  c(2, 2, 0, 0, 1, 1)
)
test_loci <- data.frame(
  name = paste0("rs", 1:6),
  chromosome = paste0("chr", c(1, 1, 1, 1, 2, 2)),
  position = as.integer(c(3, 5, 65, 343, 23, 456)),
  genetic_dist = as.double(rep(0, 6)),
  allele_ref = c("A", "T", "C", "G", "C", "T"),
  allele_alt = c("T", "C", NA, "C", "G", "A")
)

test_gt <- gen_tibble(
  x = test_genotypes,
  loci = test_loci,
  indiv_meta = test_indiv_meta,
  valid_alleles = c("A", "T", "C", "G"),
  quiet = TRUE
)

# Merge the datasets using rbind
merged_gt <- rbind(ref = example_gt, target = test_gt, flip_strand = TRUE)

merged_gt

Generate a report of what would happen to each SNP in a merge

Description

This function provides an overview of the fate of each SNP in two gen_tibble objects in the case of a merge. Only SNPs found in both objects will be kept. One object is used as a reference, and SNPs in the other dataset will be flipped and/or alleles swapped as needed. SNPs that have different alleles in the two datasets will also be dropped.

Usage

rbind_dry_run(
  ref,
  target,
  use_position = FALSE,
  flip_strand = FALSE,
  quiet = FALSE
)

Arguments

Details

Loci are matched by name, unless use_position = TRUE. If the loci names of either dataset are integers (as created by gen_tibble(names_as_int = TRUE)), loci are always matched by chromosome and position (see use_position).

Value

a list with two data.frames, named target and ref. Each data.frame has nrow() equal to the number of loci in the respective dataset, and columns id (the index of the locus in the dataset), new_id (the index of the locus in the merged dataset, NA if the locus is dropped, either because it is not found in both datasets, because its alleles do not match, or because it is ambiguous), name (the locus name, also when loci are matched by position), to_flip (boolean, loci whose strand is flipped to match the reference, only in the target data.frame), to_swap (boolean, loci whose alleles, and thus genotypes, are swapped in the merge: in target, to match the order of alleles in the reference; in ref, loci whose allele_ref is missing in both datasets, which are swapped so that the observed allele becomes allele_ref, only when flip_strand = TRUE), missing_allele and missing_allele_col (for loci with one unobserved allele (NA), the allele inferred from the other dataset, and the column it fills, allele_alt or allele_ref; both are NA if nothing was inferred), and ambiguous (boolean, strand-ambiguous loci, i.e. A/T and C/G).

Examples

example_gt <- load_example_gt("gen_tbl")

# Create a second gen_tibble to merge
test_indiv_meta <- data.frame(
  id = c("x", "y", "z"),
  population = c("pop1", "pop1", "pop2")
)
test_genotypes <- rbind(
  c(1, 1, 2, 1, 1),
  c(2, 1, 2, 0, 0),
  c(2, 2, 2, 0, 1)
)
test_loci <- data.frame(
  name = paste0("rs", 1:5),
  chromosome = paste0("chr", c(1, 1, 1, 1, 2)),
  position = as.integer(c(3, 5, 65, 343, 23)),
  genetic_dist = as.double(rep(0, 5)),
  allele_ref = c("A", "T", "C", "G", "C"),
  allele_alt = c("T", "C", NA, "C", "G")
)

test_gt <- gen_tibble(
  x = test_genotypes,
  loci = test_loci,
  indiv_meta = test_indiv_meta,
  valid_alleles = c("A", "T", "C", "G"),
  quiet = TRUE
)

# Create an rbind report using rbind_dry_run
rbind_dry_run(example_gt, test_gt, flip_strand = TRUE)

Read and structure .Q files or existing matrices as q_matrix or gt_admix objects.

Description

This function reads .Q matrix files generated by external clustering algorithms (such as ADMIXTURE) and transforms them into gt_admix objects.

Usage

read_q_files(x)

Arguments

Value

• a gt_admix object containing a list of Q matrices and a list of indices for each Q matrix separated by K

Examples

q_files_path <- system.file("extdata", "anolis", package = "tidypopgen")

admix_obj <- read_q_files(q_files_path)
summary(admix_obj)

Objects exported from other packages

Description

These objects are imported from other packages. Follow the links below to see their documentation.

generics

augment(), tidy()

ggplot2

autoplot()

Get or set row names for matrix_int_names

Description

Returns integer names if present, otherwise character dimnames. Use instead of rownames() for guaranteed dispatch.

Usage

row_names(x)

row_names(x) <- value

## Default S3 replacement method:
row_names(x) <- value

## S3 replacement method for class 'matrix_int_names'
row_names(x) <- value

Arguments

Value

Integer or character vector of row names, or NULL

The modified matrix_int_names object

See Also

Other matrix_int_names_functions: [.matrix_int_names(), as.data.frame.matrix_int_names(), col_names(), matrix_int_names(), print.matrix_int_names()

Examples

my_mat <- matrix_int_names(matrix(1:6, nrow = 2),
  row_names = c(10L, 20L),
  col_names = c(100L, 200L, 300L)
)
row_names(my_mat)
row_names(my_mat) <- c(1L, 2L)
row_names(my_mat)

Scale constructor using the distruct colours

Description

A discrete fill scale using the distruct colours from distruct_colours, in order. With more groups than colours (56), the colours are recycled, with a warning, so some groups share a colour.

Usage

scale_fill_distruct(guide = "none", ...)

Arguments

Value

a scale constructor to be used with ggplot

See Also

ggplot2::scale_fill_manual(), which this function used to wrap (it can not recycle colours).

Examples

library(ggplot2)
# Create a gen_tibble of lobster genotypes
bed_file <-
  system.file("extdata", "lobster", "lobster.bed", package = "tidypopgen")
lobsters <- gen_tibble(bed_file,
  backingfile = tempfile("lobsters"),
  quiet = TRUE
)

# Fill bars of individuals per population with the distruct colours
ggplot(lobsters, aes(x = population, fill = population)) +
  geom_bar() +
  scale_fill_distruct()

The select verb for loci

Description

An equivalent to dplyr::select() that works on the genotype column of a gen_tibble, using the mini-grammar available for tidyselect. The select-like evaluation only has access to the names of the loci (i.e. it can select only based on names, not summary statistics of those loci; look at select_loci_if() for that feature.

Usage

select_loci(.data, .sel_arg)

Arguments