Package 'tidypopgen'

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

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 gen_tibble with additional columns added to the loci tibble (accessible with show_loci(). If data is missing, a tibble of the information, with a column .rownames giving the loci names.

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

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.

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

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

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 pc)

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

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

Description

For gt_pcadapt, the following types of plots are available:

• qq: a qunatile-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

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

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

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

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

Usage

## S3 method for class 'qc_report_indiv'
autoplot(
  object,
  type = c("scatter", "relatedness"),
  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

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 freqency

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

• 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

Description

Combine method for gt_admix objects

Usage

## S3 method for class 'gt_admix'
c(...)

Arguments

Value

A gt_admix object with the combined data

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

Description

Colours in the palette used by distruct

Usage

distruct_colours

Format

A vector of 60 hex colours

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

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

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,
  quiet = FALSE
)

## S3 method for class 'character'
gen_tibble(
  x,
  ...,
  parser = c("vcfR", "cpp"),
  n_cores = 1,
  chunk_size = NULL,
  valid_alleles = c("A", "T", "C", "G"),
  missing_alleles = c("0", "."),
  backingfile = NULL,
  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,
  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. 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 counted as 0 and alternative alleles will be counted as 1.

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

Value

an object of the class gen_tbl.

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

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

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

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 = "none"
)

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

• 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)

Description

This function converts a gen_tibble to a genind object from adegenet

Usage

gt_as_genind(x)

Arguments

Value

a genind object

Description

This function converts a gen_tibble to a genlight object from adegenet

Usage

gt_as_genlight(x)

Arguments

Value

a genlight object

Description

This function writes a .geno file fom 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

NOTE that we currently read all the data into memory to write the file, so this function is not suitable for very large datasets.

Value

the path of the .geno file

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)

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. NOTE that writing to bed has been optimised for speed, but writing to ped or raw is slower, especially for large datasets.

Value

the path of the saved file

Description

This function write a VCF from a gen_tibble.

Usage

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

Arguments

Value

the path of the .vcf file

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

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 fulfils 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 unstability. 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

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 availabl plots. There are also gt_dapc_tidiers for manipulating the results. For the moment, his 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.

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 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

Description

This function prepares data for various ADMIXTOOLS 2 functions fro 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().

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 = TRUE,
  fst = TRUE,
  afprod = TRUE,
  poly_only = c("f2"),
  apply_corr = TRUE,
  n_cores = 1,
  quiet = FALSE
)

Arguments

Value

SNP metadata (invisibly)

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

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

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

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()

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.

Usage

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

Arguments

Value

A gen_tibble

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_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_randomSVD is more appropriate. Finally, there is a method specifically designed for dealing with LD in large datasets, gt_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).

Details

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

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.

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.

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

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.

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

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_frobenious, 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.

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

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.

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.

Description

This function uses a sliding-window approach to look for runs of homozygosity (or heterozygosity) in a diploid genome. This function uses the package selectRUNS, which implements an approach equivalent to the one in PLINK.

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

This function returns a data frame with all runs detected in the dataset. This data frame can then be written out to a csv file. The data frame is, in turn, the input for other functions of the detectRUNS package that create plots and produce statistics from the results (see plots and statistics functions in this manual, and/or refer to the detectRUNS vignette).

If the gen_tibble is grouped, then the grouping variable is used to fill in the group table. 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. The returned dataframe contains the following seven 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)

Examples

# run the example only if we have the package installed

if (requireNamespace("detectRUNS", quietly = TRUE)) {
  sheep_ped <- system.file("extdata", "Kijas2016_Sheep_subset.ped",
    package = "detectRUNS"
  )
  sheep_gt <- tidypopgen::gen_tibble(sheep_ped,
    backingfile = tempfile(),
    quiet = TRUE
  )
  sheep_gt <- sheep_gt %>% group_by(population)
  sheep_roh <- gt_roh_window(sheep_gt)
  detectRUNS::plot_Runs(runs = sheep_roh)
}

Description

Save the tibble (and update the backing files). The gen_tibble object is saved to a file with extension .gt, togethe 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()

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

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 the function snmf from R package LEA.

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)

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. Tests for this function are in test_gt_order_loci.R

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().

Value

a gen_tibble with a backing file (i.e. a new File Backed Matrix)

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

Description

Estimate observed heterozygosity (H_obs) for each individual (i.e. the frequency of loci that are heterozygous in an individual).

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

a vector of heterozygosities, one per individuals in the gen_tibble

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_df'
indiv_inbreeding(.x, method = c("WG17"), allele_sharing_mat = NULL, ...)

Arguments

Value

a numeric vector of inbreeding coefficients.

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 heterozygosities, one per individuals in the gen_tibble

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

Description

This functions checks that all SNPs in a chromosome are adjacent in the loci table, and that positions are sorted within chromosomes.

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

a logical vector defining which loci are transversions

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.

Usage

loci_alt_freq(.x, n_cores, block_size, ...)

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

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

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

loci_maf(.x, n_cores, block_size, ...)

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

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

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

Arguments

Value

a vector of frequencies, one per locus

Description

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

Usage

loci_chromosomes(.x, ...)

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

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

Arguments

Value

a character vector of chromosomes

Description

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

Usage

loci_hwe(.x, ...)

## S3 method for class 'tbl_df'
loci_hwe(.x, mid_p = TRUE, ...)

## S3 method for class 'vctrs_bigSNP'
loci_hwe(.x, mid_p = TRUE, ...)

## S3 method for class 'grouped_df'
loci_hwe(.x, ...)

Arguments

Details

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

NOTE There are no tests for this function yet! Unit tests are needed.

Value

a vector of probabilities from HWE exact test, one per locus

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).

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.

Usage

loci_ld_clump(.x, ...)

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

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

## S3 method for class 'grouped_df'
loci_ld_clump(.x, ...)

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')

Description

Estimate the rate of missingness at each locus.

Usage

loci_missingness(
  .x,
  as_counts = FALSE,
  n_cores = bigstatsr::nb_cores(),
  block_size,
  ...
)

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

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

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

Arguments

Value

a vector of frequencies, one per locus

Description

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

Usage

loci_names(.x, ...)

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

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

Arguments

Value

a character vector of names

Description

Use the loci table to define which loci are transitions

Usage

loci_transitions(.x, ...)

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

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

## S3 method for class 'grouped_df'
loci_transitions(.x, ...)

Arguments

Value

a logical vector defining which loci are transitions

Description

Use the loci table to define which loci are transversions

Usage

loci_transversions(.x, ...)

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

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

## S3 method for class 'grouped_df'
loci_transversions(.x, ...)

Arguments

Value

a logical vector defining which loci are transversions

Description

This function computes the Allele Sharing matrix. Estimates Allele Sharing (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))
)

Arguments

Value

a matrix of allele sharing between all pairs of individuals

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 follwng 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))
)

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

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))
)

Arguments

Details

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

Value

a bigstatsr::FBM of proportion or adjusted counts, or a list of two bigstatsr::FBM matrices, one of counts of IBS by alleles, and one of number of valid alleles (i.e. 2n_loci - 2missing_loci)

Description

This function computes the KING-robust estimator of kinship.

Usage

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

Arguments

Details

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

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.

Description

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

• 'Hudson': Hudson's formulation, as derived in Bhatia et al (2013) for diploids.

• 'Nei87' : Fst according to Nei (1987) - includes the correction for heterozygosity when computing Ht, and is equivalent to hierfstat::pairwise.neifst(),

• 'WC84' : Weir and Cockerham (1984), correcting for missing data and is equivalent to hierfstat::pairwise.WCfst().

Usage

pairwise_pop_fst(
  .x,
  tidy = TRUE,
  by_locus = FALSE,
  method = c("Hudson", "Nei87", "WC84"),
  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.

Value

if tidy=TRUE, 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 tidy=FALSE, 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.

Description

This function computes population specific FIS, using either the approach of Nei 1987 (as computed by hierfstat::basic.stats()) or of Weir and Goudet 2017 (as computed by hierfstat::fis.dosage()).

Usage

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

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

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

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 $k$th 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 monorphic 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.

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

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 $k$th 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.

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

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,
  ...
)

Arguments

Value

a matrix of predictions, 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.

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

Description

Return QC information to assess loci (Observed heterozygosity and missingness).

Usage

qc_report_indiv(.x, ...)

## S3 method for class 'tbl_df'
qc_report_indiv(.x, kings_threshold = NULL, ...)

## S3 method for class 'grouped_df'
qc_report_indiv(.x, kings_threshold = NULL, ...)

Arguments

Value

a tibble with 2 elements: het_obs and missingness

Description

Return QC information to assess loci (MAF, missingness and HWE test).

Usage

qc_report_loci(.x, ...)

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

## S3 method for class 'grouped_df'
qc_report_loci(.x, ...)

Arguments

Value

a tibble with 3 elements: maf, missingness and hwe_p

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

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, a column id with the locus name, and boolean columns to_keep (the valid loci that will be kept in the merge), alleles_mismatched (loci found in both datasets but with mismatched alleles, leading to those loci being dropped), to_flip (loci that need to be flipped to align the two datasets, only found in target data.frame) and to_swap (loci for which the order of alleles needs to be swapped to align the two datasets, target data.frame)

Description

This function combined two gen_tibbles. By defaults, it subsets the loci and swaps ref and alt alleles to make the two datasets compatible (this behaviour can be switched off with as_is). 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

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.

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

Description

A wrapper around ggplot2::scale_fill_manual(), using the distruct colours from distruct_colours.

Usage

scale_fill_distruct(guide = "none", ...)

Arguments

Value

a scale constructor to be used with ggplot

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

Details

Note that the select_loci verb does not modify the backing FBM files, but rather it subsets the list of loci to be used stored in the gen_tibble.

Value

a gen_tibble with a subset of the loci.

Description

An equivalent to dplyr::select_if() that works on the genotype column of a gen_tibble. This function has access to the genotypes (and thus can work on summary statistics to select), but not the names of the loci (look at select_loci() for that feature.

Usage

select_loci_if(.data, .sel_logical)

Arguments

Details

#' Note that the select_loci_if verb does not modify the backing FBM files, but rather it subsets the list of loci to be used stored in the gen_tibble.

Value

a subset of the list of loci in the gen_tibble

Description

Extract the genotypes (as a matrix) from a gen_tibble.

Usage

show_genotypes(.x, indiv_indices = NULL, loci_indices = NULL, ...)

## S3 method for class 'tbl_df'
show_genotypes(.x, indiv_indices = NULL, loci_indices = NULL, ...)

## S3 method for class 'vctrs_bigSNP'
show_genotypes(.x, indiv_indices = NULL, loci_indices = NULL, ...)

Arguments

Value

a matrix of counts of the alternative alleles (see show_loci()) to extract information on the alleles for those loci from a gen_tibble.

Description

Extract and set the information on loci from a gen_tibble.

Usage

show_loci(.x, ...)

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

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

show_loci(.x) <- value

## S3 replacement method for class 'tbl_df'
show_loci(.x) <- value

## S3 replacement method for class 'vctrs_bigSNP'
show_loci(.x) <- value

Arguments

Value

a tibble::tibble of information (see gen_tibble for details on compulsory columns that will always be present)

Description

Extract the ploidy information from a gen_tibble. NOTE that this function does not return the ploidy level for each individual (that is obtained with indiv_ploidy); instead, it returns an integer which is either the ploidy level of all individuals (e.g. 2 indicates all individuals are diploid), or a 0 to indicate mixed ploidy.

Usage

show_ploidy(.x, ...)

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

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

Arguments

Value

the ploidy (0 indicates mixed ploidy)

Description

This function computes the Allele Sharing matrix. Estimates Allele Sharing (matching in hierfstat)) 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

snp_allele_sharing(
  X,
  ind.row = bigstatsr::rows_along(X),
  ind.col = bigstatsr::cols_along(X),
  block.size = bigstatsr::block_size(nrow(X))
)

Arguments

Value

a matrix of allele sharing between all pairs of individuals

Description

This function computes the IBS matrix.

Usage

snp_ibs(
  X,
  ind.row = bigstatsr::rows_along(X),
  ind.col = bigstatsr::cols_along(X),
  type = c("proportion", "adjusted_counts", "raw_counts"),
  block.size = bigstatsr::block_size(nrow(X))
)

Arguments

Details

Note that monomorphic sites are currently counted. Should we filter them beforehand? What does plink do?

Value

if as.counts = TRUE function returns a list of two bigstatsr::FBM matrices, one of counts of IBS by alleles (i.e. 2*n loci), and one of valid alleles (i.e. 2 * n_loci - 2 * missing_loci). If as.counts = FALSE returns a single matrix of IBS proportions.

Description

This function computes the KING-robust estimator of kinship.

Usage

snp_king(
  X,
  ind.row = bigstatsr::rows_along(X),
  ind.col = bigstatsr::cols_along(X),
  block.size = bigstatsr::block_size(nrow(X)) * 4
)

Arguments

Value

a square symmetrical matrix of relationship coefficients between individuals

Description

Summary method for gt_admix objects

Usage

## S3 method for class 'gt_admix'
summary(object, ...)

Arguments

Value

A summary of the gt_admix object

Description

This function creates a summary of the merge report generated by rbind_dry_run()

Usage

## S3 method for class 'rbind_report'
summary(object, ..., ref_label = "reference", target_label = "target")

Arguments

Value

NULL (prints a summary to the console)

Description

A theme to remove most plot decorations, matching the look of plots created with distruct.

Usage

theme_distruct()

Value

a ggplot2::theme

Description

This summarizes information about the components of a gt_dapc from the tidypopgen package. The parameter matrix determines which element is returned.

Usage

## S3 method for class 'gt_dapc'
tidy(x, matrix = "eigenvalues", ...)

Arguments

Value

A tibble::tibble with columns depending on the component of DAPC being tidied.

If "scores" each row in the tidied output corresponds to the original data in PCA space. The columns are:

If matrix is "loadings", each row in the tidied output corresponds to information about the principle components in the original space. The columns are:

If "eigenvalues", the columns are:

See Also

gt_dapc() augment.gt_dapc()

Description

This summarizes information about the components of a gt_pca from the tidypopgen package. The parameter matrix determines which element is returned. Column names of the tidied output match those returned by broom::tidy.prcomp, the tidier for the standard PCA objects returned by stats::prcomp.

Usage

## S3 method for class 'gt_pca'
tidy(x, matrix = "eigenvalues", ...)

Arguments

Value

A tibble::tibble with columns depending on the component of PCA being tidied.

If "scores" each row in the tidied output corresponds to the original data in PCA space. The columns are:

If matrix is "loadings", each row in the tidied output corresponds to information about the principle components in the original space. The columns are:

If "eigenvalues", the columns are:

See Also

gt_pca_autoSVD() augment_gt_pca

Description

Takes a q_matrix object, which is a matrix, and returns a tidied tibble.

Usage

## S3 method for class 'q_matrix'
tidy(x, data, ...)

Arguments

Value

A tidied tibble containing columns: