From 3be63082259a9322aa23668f4c604da355c81e83 Mon Sep 17 00:00:00 2001 From: Andrea Manica Date: Fri, 9 Aug 2024 10:45:46 +0100 Subject: [PATCH] Fix test --- NAMESPACE | 1 + R/gt_as_vcf.R | 4 ++- R/gt_pcadapt.R | 48 ++++++++++++++++++++++++++++++++ man/gt_pcadapt.Rd | 30 ++++++++++++++++++++ tests/testthat/test_gt_pcadapt.R | 10 +++++++ 5 files changed, 92 insertions(+), 1 deletion(-) create mode 100644 R/gt_pcadapt.R create mode 100644 man/gt_pcadapt.Rd create mode 100644 tests/testthat/test_gt_pcadapt.R diff --git a/NAMESPACE b/NAMESPACE index b918aed0..f758402c 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -100,6 +100,7 @@ export(gt_load) export(gt_pca_autoSVD) export(gt_pca_partialSVD) export(gt_pca_randomSVD) +export(gt_pcadapt) export(gt_roh_window) export(gt_save) export(gt_set_imputed) diff --git a/R/gt_as_vcf.R b/R/gt_as_vcf.R index 4ed65713..b7178feb 100644 --- a/R/gt_as_vcf.R +++ b/R/gt_as_vcf.R @@ -13,9 +13,11 @@ #' gt_as_vcf <- function(x, file = NULL, chunk_size = NULL, overwrite = FALSE){ + # check that x is a gen_tibble + stopifnot_gen_tibble(x) # vcf writing currently only works for diploid data - stopifnot_diploid(x) + stopifnot_diploid(x$genotypes) if (is.null(file)){ file <- bigstatsr::sub_bk(attr(x$genotypes,"bigsnp")$genotypes$backingfile,paste0(".vcf")) diff --git a/R/gt_pcadapt.R b/R/gt_pcadapt.R new file mode 100644 index 00000000..6c5c6932 --- /dev/null +++ b/R/gt_pcadapt.R @@ -0,0 +1,48 @@ +#' pcadapt analysis on a `gen_tibble` object +#' +#' 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)](https://doi.org/10.1534/genetics.116.195214), +#' See the R package `pcadapt`, which provides extensive +#' documentation and examples. +#' +#' Internally, this function uses the `snp_pcadapt` function from the `bigsnpr` package. +#' @param x A `gen_tibble` object. +#' @param pca a [`gt_pca`] object, as returned by `gt_pca_partialSVD()` or `gt_pca_randomSVD()`. +#' @param k Number of principal components to use in the analysis. +#' @param n_cores Number of cores to use. +#' @returns An object of subclass `gt_pcadapt`, a subclass of `mhtest`. +#' @export + +gt_pcadapt <- function(x, pca, k, n_cores = 1) { + stopifnot_gen_tibble(x) + if (!inherits(pca, "gt_pca")) { + stop("pca must be a gt_pca object") + } + # check that k is a scalar + if (!is.numeric(k) || length(k) != 1) { + stop("k must be a scalar") + } + # check that k is not larger than the number of components in pca + if (k > ncol(pca$v)) { + stop("K must be less than or equal to the number of components in pca") + } + # set imputation if needed + if (gt_has_imputed(x) && gt_uses_imputed(x)==FALSE){ + gt_set_imputed(x, set = TRUE) + on.exit(gt_set_imputed(x, set = FALSE)) + } + + # Run the analysis + res <- bigsnpr::snp_pcadapt( + G = .gt_get_bigsnp(x)$genotypes, + U.row = pca$u[, 1:k, drop = FALSE], + ind.row = .gt_bigsnp_rows(x), + ind.col = .gt_bigsnp_cols(x), + ncores = n_cores + ) + + class(res) <- c("gt_pcadapt", class(res)) + # Return the result + return(res) +} diff --git a/man/gt_pcadapt.Rd b/man/gt_pcadapt.Rd new file mode 100644 index 00000000..a56421e4 --- /dev/null +++ b/man/gt_pcadapt.Rd @@ -0,0 +1,30 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/gt_pcadapt.R +\name{gt_pcadapt} +\alias{gt_pcadapt} +\title{pcadapt analysis on a \code{gen_tibble} object} +\usage{ +gt_pcadapt(x, pca, k, n_cores = 1) +} +\arguments{ +\item{x}{A \code{gen_tibble} object.} + +\item{pca}{a \code{\link{gt_pca}} object, as returned by \code{gt_pca_partialSVD()} or \code{gt_pca_randomSVD()}.} + +\item{k}{Number of principal components to use in the analysis.} + +\item{n_cores}{Number of cores to use.} +} +\value{ +An object of subclass \code{gt_pcadapt}, a subclass of \code{mhtest}. +} +\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 \href{https://doi.org/10.1534/genetics.116.195214}{Luu et al. (2017)}, +See the R package \code{pcadapt}, which provides extensive +documentation and examples. +} +\details{ +Internally, this function uses the \code{snp_pcadapt} function from the \code{bigsnpr} package. +} diff --git a/tests/testthat/test_gt_pcadapt.R b/tests/testthat/test_gt_pcadapt.R new file mode 100644 index 00000000..ce726b9a --- /dev/null +++ b/tests/testthat/test_gt_pcadapt.R @@ -0,0 +1,10 @@ +test_that("gt_pcadapt works on gt_pca object",{ + bed_file <- system.file("extdata", "example-missing.bed", package = "bigsnpr") + missing_gt <- gen_tibble(bed_file, backingfile = tempfile("missing_"),quiet=TRUE) + missing_gt <- gt_impute_simple(missing_gt, method = "mode") + missing_pca <- missing_gt %>% gt_pca_partialSVD() + missing_pcadapt <- gt_pcadapt(missing_gt, missing_pca, k = 3) + expect_true((inherits(missing_pcadapt, "gt_pcadapt"))) +}) + +