diff --git a/DESCRIPTION b/DESCRIPTION index 6f5bfc3..cfa805a 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: AllelicSeries Title: Allelic Series Test -Version: 0.0.4.1 +Version: 0.0.5.0 Authors@R: c(person(given = "Zachary", family = "McCaw", @@ -14,7 +14,7 @@ Authors@R: comment = c(ORCID = "0000-0002-1748-625X")), person(given = "insitro", role = c("cph")) ) -Description: Implementation of gene-level rare variant association tests targeting allelic series: genes where increasingly deleterious mutations have increasingly large phenotypic effects. The COding-variant Allelic Series Test (COAST) operates on the benign missense variants (BMVs), deleterious missense variants (DMVs), and protein truncating variants (PTVs) within a gene. COAST uses a set of adjustable weights that tailor the test towards rejecting the null hypothesis for genes where the average magnitude of effect increases monotonically from BMVs to DMVs to PTVs. See McCaw ZR, O’Dushlaine C, Somineni H, Bereket M, Klein C, Karaletsos T, Casale FP, Koller D, Soare TW. (2022) "An allelic series rare variant association test for candidate gene discovery" . +Description: Implementation of gene-level rare variant association tests targeting allelic series: genes where increasingly deleterious mutations have increasingly large phenotypic effects. The COding-variant Allelic Series Test (COAST) operates on the benign missense variants (BMVs), deleterious missense variants (DMVs), and protein truncating variants (PTVs) within a gene. COAST uses a set of adjustable weights that tailor the test towards rejecting the null hypothesis for genes where the average magnitude of effect increases monotonically from BMVs to DMVs to PTVs. See McCaw ZR, O’Dushlaine C, Somineni H, Bereket M, Klein C, Karaletsos T, Casale FP, Koller D, Soare TW. (2023) "An allelic series rare variant association test for candidate gene discovery" . License: BSD_3_clause + file LICENSE Encoding: UTF-8 Imports: @@ -25,7 +25,7 @@ LinkingTo: Rcpp, RcppArmadillo Roxygen: list(markdown = TRUE) -RoxygenNote: 7.2.0 +RoxygenNote: 7.2.3 Suggests: knitr, rmarkdown, diff --git a/NAMESPACE b/NAMESPACE index 6575cab..018120d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -5,6 +5,7 @@ export(ASKAT) export(Aggregator) export(COAST) export(Comparator) +export(CountAlleles) export(DGP) export(OLS) importFrom(Rcpp,sourceCpp) diff --git a/NEWS.md b/NEWS.md new file mode 100644 index 0000000..e0cd170 --- /dev/null +++ b/NEWS.md @@ -0,0 +1,6 @@ +## Version 0.5.0 + +* Added option (`min_mac`) to filter the variant set to only include those variants having at least a minimum minor allele count (10 is recommended). +* Added a function (`CountAlleles`) to count the number of alleles of each variant category present in the genotype matrix. Also allows for counting the number of carriers of each type of allele. +* By default, `COAST` now reports the number of alleles of each variant category that contributed to the test. + diff --git a/R/allelic_series_test.R b/R/allelic_series_test.R index 4764194..95c0939 100644 --- a/R/allelic_series_test.R +++ b/R/allelic_series_test.R @@ -1,5 +1,5 @@ # Purpose: Allelic series test. -# Updated: 2023-10-04 +# Updated: 2024-04-03 # Default weights. DEFAULT_WEIGHTS <- c(1, 2, 3) @@ -14,6 +14,7 @@ DEFAULT_WEIGHTS <- c(1, 2, 3) #' @param indicator Convert raw counts to indicators? Default: FALSE. #' @param method Method for aggregating across categories: #' {"none", "max", "sum"}. Default: "none". +#' @param min_mac Minimum minor allele count for inclusion. Default: 0. #' @param weights Annotation category weights. #' @return (n x 3) Numeric matrix without weighting, (n x 1) numeric matrix #' with weighting. @@ -24,8 +25,18 @@ Aggregator <- function( drop_empty = TRUE, indicator = FALSE, method = "none", + min_mac = 0, weights = DEFAULT_WEIGHTS ) { + + # Minor allele count filtering. + if (min_mac > 0) { + mac <- apply(geno, 2, sum) + keep <- (mac >= min_mac) + + anno <- anno[keep] + geno <- geno[, keep, drop = FALSE] + } # Sum to categories. bmv <- apply(geno[, anno == 0, drop = FALSE], 1, sum) @@ -95,6 +106,7 @@ Aggregator <- function( #' @param is_pheno_binary Is the phenotype binary? Default: FALSE. #' @param method Method for aggregating across categories: {"none", "max", #' "sum"}. Default: "none". +#' @param min_mac Minimum minor allele count for inclusion. Default: 0. #' @param score_test Run a score test? If FALSE, performs a Wald test. #' @param weights (3 x 1) annotation category weights. #' @return Numeric p-value. @@ -120,6 +132,7 @@ ASBT <- function( indicator = FALSE, is_pheno_binary = FALSE, method = "none", + min_mac = 0, score_test = FALSE, weights = DEFAULT_WEIGHTS ) { @@ -152,6 +165,7 @@ ASBT <- function( drop_empty = TRUE, indicator = indicator, method = method, + min_mac = min_mac, weights = weights ) @@ -191,6 +205,7 @@ ASBT <- function( #' Default: TRUE. Ignored if phenotype is binary. #' @param covar (n x p) covariate matrix. Defaults to an (n x 1) intercept. #' @param is_pheno_binary Is the phenotype binary? Default: FALSE. +#' @param min_mac Minimum minor allele count for inclusion. Default: 0. #' @param return_null_model Return the null model in addition to the p-value? #' Useful if running additional SKAT tests. Default: FALSE. #' @param weights (3 x 1) annotation category weights. @@ -216,6 +231,7 @@ ASKAT <- function( apply_int = TRUE, covar = NULL, is_pheno_binary = FALSE, + min_mac = 0, return_null_model = FALSE, weights = DEFAULT_WEIGHTS ) { @@ -241,6 +257,15 @@ ASKAT <- function( weights = weights ) + # Minor allele count filtering. + if (min_mac >= 0) { + mac <- apply(geno, 2, sum) + keep <- (mac > min_mac) + + anno <- anno[keep] + geno <- geno[, keep, drop = FALSE] + } + # Alternate allele frequencies. aaf <- apply(geno, 2, mean) / 2 @@ -320,6 +345,8 @@ ASKAT <- function( #' @param include_orig_skato_ptv Include the original version of SKAT-O applied #' to PTV variants only in the omnibus test? Default: FALSE. #' @param is_pheno_binary Is the phenotype binary? Default: FALSE. +#' @param min_mac Minimum minor allele count for inclusion. Default: 0. +#' @param return_counts Include minor allele counts in output? Default: TRUE. #' @param return_omni_only Return only the omnibus p-value? Default: FALSE. #' @param score_test Use a score test for burden analysis? If FALSE, uses a #' Wald test. @@ -348,6 +375,8 @@ COAST <- function( include_orig_skato_all = FALSE, include_orig_skato_ptv = FALSE, is_pheno_binary = FALSE, + min_mac = 0, + return_counts = TRUE, return_omni_only = FALSE, score_test = FALSE, weights = DEFAULT_WEIGHTS @@ -381,6 +410,7 @@ COAST <- function( geno = geno, pheno = pheno, apply_int = apply_int, + min_mac = min_mac, is_pheno_binary = is_pheno_binary, score_test = score_test, ... @@ -427,6 +457,7 @@ COAST <- function( geno = geno, pheno = pheno, is_pheno_binary = is_pheno_binary, + min_mac = min_mac, return_null_model = TRUE, weights = weights ) @@ -463,5 +494,11 @@ COAST <- function( } else { out <- c(p_val, p_omni = p_omni) } + + if (return_counts) { + counts <- CountAlleles(anno = anno, geno = geno, min_mac = min_mac) + out <- c(counts, out) + } + return(out) } diff --git a/R/utilities.R b/R/utilities.R index 64c7300..dadb60f 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -1,5 +1,54 @@ # Purpose: Utility functions. -# Updated: 2023-04-25 +# Updated: 2024-04-03 + +#' Count Alleles +#' +#' Count the number of non-zero alleles bearing each variant annotation. +#' +#' @param anno (snps x 1) annotation vector with values in c(0, 1, 2). +#' @param geno (n x snps) genotype matrix. +#' @param count_carriers If true, counts the number of carriers rather than +#' the number of alleles. +#' @param min_mac Minimum minor allele count for inclusion. Default: 0. +#' @return 3 x 1 numeric vector with the counts of BMVs, DMVs, and PTVs. +#' @export +CountAlleles <- function( + anno, + geno, + count_carriers = FALSE, + min_mac = 0 +) { + + # Minor allele count filtering. + if (min_mac > 0) { + mac <- apply(geno, 2, sum) + keep <- (mac >= min_mac) + + anno <- anno[keep] + geno <- geno[, keep, drop = FALSE] + } + + # Category counts. + if (count_carriers) { + + # Count carriers. + n_bmv <- sum(apply(geno[, anno == 0, drop = FALSE], 1, sum) > 0) + n_dmv <- sum(apply(geno[, anno == 1, drop = FALSE], 1, sum) > 0) + n_ptv <- sum(apply(geno[, anno == 2, drop = FALSE], 1, sum) > 0) + + } else { + + # Count alleles. + n_bmv <- sum(geno[, anno == 0]) + n_dmv <- sum(geno[, anno == 1]) + n_ptv <- sum(geno[, anno == 2]) + + } + + # Output. + out <- c(n_bmv = n_bmv, n_dmv = n_dmv, n_ptv = n_ptv) + return(out) +} #' Linear Association Test diff --git a/README.md b/README.md index 2d83a38..4026ff4 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,20 @@ # Allelic Series -Implementation of gene-level rare variant association tests targeting allelic series: genes where increasingly deleterious mutations have increasingly large phenotypic effects. The COding-variant Allelic Series Test (COAST) operates on the benign missense variants (BMVs), deleterious missense variants (DMVs), and protein truncating variants (PTVs) within a gene. COAST uses a set of adjustable weights that tailor the test towards rejecting the null hypothesis for genes where the average magnitude of effect increases monotonically from BMVs to DMVs to PTVs. See McCaw ZR, O’Dushlaine C, Somineni H, Bereket M, Klein C, Karaletsos T, Casale FP, Koller D, Soare TW. (2023) "An allelic-series rare-variant association test for candidate-gene discovery" [doi:10.1016/j.ajhg.2023.07.001](https://www.cell.com/ajhg/fulltext/S0002-9297(23)00241-0). +This package implements gene-level rare variant association tests +targeting allelic series: genes where increasingly deleterious mutations +have increasingly large phenotypic effects. The main COding-variant +Allelic Series Test (COAST) operates on the benign missense variants +(BMVs), deleterious missense variants (DMVs), and protein truncating +variants (PTVs) within a gene. COAST uses a set of adjustable weights +that tailor the test towards rejecting the null for genes where the +average magnitude of phenotypic effect increases monotonically from BMVs +to DMVs to PTVs. Such genes are of candidate therapeutic interest due to +the existence of a dose-response relationship between gene functionality +and phenotypic impact. See McCaw ZR, O’Dushlaine C, Somineni H, Bereket +M, Klein C, Karaletsos T, Casale FP, Koller D, Soare TW. (2023) “An +allelic-series rare-variant association test for candidate-gene +discovery” +[doi:10.1016/j.ajhg.2023.07.001](https://www.cell.com/ajhg/fulltext/S0002-9297(23)00241-0). # Installation @@ -28,11 +42,13 @@ The example `data` are a list with the following components: - `anno`: An `snps` by 1 annotation vector coded as 0 for benign missense variants (BMVs), 1 for deleterious missense variants (DMVs), - and 2 for protein truncating variants (PTVs). + and 2 for protein truncating variants (PTVs). Note that the values of + (0, 1, 2) simply identify different categories of variants; `weights` + other than these can be set when performing the association test. - `covar`: An `n` by 6 covariate matrix including an intercept `int`, - and covariates representing `age`, `sex`, and 3 genetic PCs (`pc1` … - `pc3`). + and covariates representing `age`, `sex`, and 3 genetic PCs (`pc1`, + `pc2`, `pc3`). - `geno`: An `n` by `snps` genotype matrix with additive coding and minor allele frequencies between 0.5% and 1.0%. @@ -67,24 +83,29 @@ genotype matrix, and the phenotype vector. included manually, if desired. - `weights` encodes the relative importance of BMVs, DMVs, and PTVs. The - example weights of `c(1, 2, 3)` target a genetic architecture where effect sizes increase with increasing deleteriousness: - BMVs have an effect of 1, DMVs have an effect of 2, and PTVs have an effect of 3. Weights of - `c(1, 1, 1)` target instead a genetic architecture where all variant - types have equivalent expected magnitudes. + example weights of `c(1, 2, 3)` target a genetic architecture where + effect sizes increase with increasing deleteriousness: BMVs have an + effect of 1, DMVs have an effect of 2, and PTVs have an effect of 3. + Weights of `c(1, 1, 1)` target instead a genetic architecture where + all variant types have equivalent expected magnitudes. ``` r show(results) ``` - ## p_count p_ind p_max_count p_max_ind p_sum_count - ## 7.707024e-29 6.745269e-06 4.299938e-13 6.228669e-07 6.756953e-18 - ## p_sum_ind p_allelic_skat p_omni - ## 1.894500e-06 8.507977e-08 9.248429e-28 + ## n_bmv n_dmv n_ptv p_count p_ind + ## 2.870000e+02 1.620000e+02 6.100000e+01 3.112702e-26 1.322084e-09 + ## p_max_count p_max_ind p_sum_count p_sum_ind p_allelic_skat + ## 3.076876e-10 5.374363e-09 1.661854e-20 2.554417e-11 2.658137e-07 + ## p_omni + ## 3.735235e-25 -By default, the output of `COAST` is a vector of p-values, corresponding -to the different components of the allelic series test and the overall -omnibus test (`p_omni`). To return the omnibus p-value only, specify -`return_omni_only = TRUE` when calling `COAST`. +By default, the output of `COAST` includes counts for the number of +alleles of each variant class that contributed to the test, and a vector +of p-values, corresponding to the different components of the allelic +series test. The final, overall p-value is given by `p_omni`. To return +the omnibus p-value only, specify `return_omni_only = TRUE` when calling +`COAST`. ## Robust omnibus test @@ -112,11 +133,17 @@ results <- COAST( show(results) ``` - ## p_count p_ind p_max_count p_max_ind p_sum_count - ## 7.707024e-29 6.745269e-06 4.299938e-13 6.228669e-07 6.756953e-18 - ## p_sum_ind p_allelic_skat p_orig_skat_all p_orig_skat_ptv p_omni - ## 1.894500e-06 8.507977e-08 1.250426e-05 2.237988e-13 9.248429e-28 + ## n_bmv n_dmv n_ptv p_count p_ind + ## 2.870000e+02 1.620000e+02 6.100000e+01 3.112702e-26 1.322084e-09 + ## p_max_count p_max_ind p_sum_count p_sum_ind p_allelic_skat + ## 3.076876e-10 5.374363e-09 1.661854e-20 2.554417e-11 2.658137e-07 + ## p_orig_skat_all p_orig_skat_ptv p_omni + ## 1.548324e-05 6.632119e-08 3.735235e-25 ## Loading genotypes -The [genio](https://CRAN.R-project.org/package=genio) and [rbgen](https://enkre.net/cgi-bin/code/bgen/wiki?name=rbgen) packages may be used to load PLINK and BGEN genotypes in R, respectively. Moreover, [PLINK](https://www.cog-genomics.org/plink/2.0/) enables conversion between the file types. +The [genio](https://CRAN.R-project.org/package=genio) and +[rbgen](https://enkre.net/cgi-bin/code/bgen/wiki?name=rbgen) packages +may be used to load PLINK and BGEN genotypes in R, respectively. +Moreover, [PLINK](https://www.cog-genomics.org/plink/2.0/) enables +conversion between these file types. diff --git a/man/ASBT.Rd b/man/ASBT.Rd index 0609e82..d034ca8 100644 --- a/man/ASBT.Rd +++ b/man/ASBT.Rd @@ -13,6 +13,7 @@ ASBT( indicator = FALSE, is_pheno_binary = FALSE, method = "none", + min_mac = 0, score_test = FALSE, weights = DEFAULT_WEIGHTS ) @@ -36,6 +37,8 @@ Default: TRUE. Ignored if phenotype is binary.} \item{method}{Method for aggregating across categories: {"none", "max", "sum"}. Default: "none".} +\item{min_mac}{Minimum minor allele count for inclusion. Default: 0.} + \item{score_test}{Run a score test? If FALSE, performs a Wald test.} \item{weights}{(3 x 1) annotation category weights.} diff --git a/man/ASKAT.Rd b/man/ASKAT.Rd index c00f722..3c3540f 100644 --- a/man/ASKAT.Rd +++ b/man/ASKAT.Rd @@ -11,6 +11,7 @@ ASKAT( apply_int = TRUE, covar = NULL, is_pheno_binary = FALSE, + min_mac = 0, return_null_model = FALSE, weights = DEFAULT_WEIGHTS ) @@ -29,6 +30,8 @@ Default: TRUE. Ignored if phenotype is binary.} \item{is_pheno_binary}{Is the phenotype binary? Default: FALSE.} +\item{min_mac}{Minimum minor allele count for inclusion. Default: 0.} + \item{return_null_model}{Return the null model in addition to the p-value? Useful if running additional SKAT tests. Default: FALSE.} diff --git a/man/Aggregator.Rd b/man/Aggregator.Rd index 7289681..cea47bd 100644 --- a/man/Aggregator.Rd +++ b/man/Aggregator.Rd @@ -10,6 +10,7 @@ Aggregator( drop_empty = TRUE, indicator = FALSE, method = "none", + min_mac = 0, weights = DEFAULT_WEIGHTS ) } @@ -25,6 +26,8 @@ Aggregator( \item{method}{Method for aggregating across categories: {"none", "max", "sum"}. Default: "none".} +\item{min_mac}{Minimum minor allele count for inclusion. Default: 0.} + \item{weights}{Annotation category weights.} } \value{ diff --git a/man/AllelicSeries-package.Rd b/man/AllelicSeries-package.Rd index 32669fc..513f2ac 100644 --- a/man/AllelicSeries-package.Rd +++ b/man/AllelicSeries-package.Rd @@ -6,7 +6,7 @@ \alias{AllelicSeries-package} \title{AllelicSeries: Allelic Series Test} \description{ -Implementation of gene-level rare variant association tests targeting allelic series: genes where increasingly deleterious mutations have increasingly large phenotypic effects. The COding-variant Allelic Series Test (COAST) operates on the benign missense variants (BMVs), deleterious missense variants (DMVs), and protein truncating variants (PTVs) within a gene. COAST uses a set of adjustable weights that tailor the test towards rejecting the null hypothesis for genes where the average magnitude of effect increases monotonically from BMVs to DMVs to PTVs. See McCaw ZR, O’Dushlaine C, Somineni H, Bereket M, Klein C, Karaletsos T, Casale FP, Koller D, Soare TW. (2022) "An allelic series rare variant association test for candidate gene discovery" \doi{10.1101/2022.12.23.521658}. +Implementation of gene-level rare variant association tests targeting allelic series: genes where increasingly deleterious mutations have increasingly large phenotypic effects. The COding-variant Allelic Series Test (COAST) operates on the benign missense variants (BMVs), deleterious missense variants (DMVs), and protein truncating variants (PTVs) within a gene. COAST uses a set of adjustable weights that tailor the test towards rejecting the null hypothesis for genes where the average magnitude of effect increases monotonically from BMVs to DMVs to PTVs. See McCaw ZR, O’Dushlaine C, Somineni H, Bereket M, Klein C, Karaletsos T, Casale FP, Koller D, Soare TW. (2023) "An allelic series rare variant association test for candidate gene discovery" \doi{10.1016/j.ajhg.2023.07.001}. } \author{ \strong{Maintainer}: Zachary McCaw \email{zmccaw@insitro.com} (\href{https://orcid.org/0000-0002-2006-9828}{ORCID}) diff --git a/man/COAST.Rd b/man/COAST.Rd index 7066c54..37dc654 100644 --- a/man/COAST.Rd +++ b/man/COAST.Rd @@ -13,6 +13,8 @@ COAST( include_orig_skato_all = FALSE, include_orig_skato_ptv = FALSE, is_pheno_binary = FALSE, + min_mac = 0, + return_counts = TRUE, return_omni_only = FALSE, score_test = FALSE, weights = DEFAULT_WEIGHTS @@ -38,6 +40,10 @@ to PTV variants only in the omnibus test? Default: FALSE.} \item{is_pheno_binary}{Is the phenotype binary? Default: FALSE.} +\item{min_mac}{Minimum minor allele count for inclusion. Default: 0.} + +\item{return_counts}{Include minor allele counts in output? Default: TRUE.} + \item{return_omni_only}{Return only the omnibus p-value? Default: FALSE.} \item{score_test}{Use a score test for burden analysis? If FALSE, uses a diff --git a/man/CountAlleles.Rd b/man/CountAlleles.Rd new file mode 100644 index 0000000..bb3169d --- /dev/null +++ b/man/CountAlleles.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utilities.R +\name{CountAlleles} +\alias{CountAlleles} +\title{Count Alleles} +\usage{ +CountAlleles(anno, geno, count_carriers = FALSE, min_mac = 0) +} +\arguments{ +\item{anno}{(snps x 1) annotation vector with values in c(0, 1, 2).} + +\item{geno}{(n x snps) genotype matrix.} + +\item{count_carriers}{If true, counts the number of carriers rather than +the number of alleles.} + +\item{min_mac}{Minimum minor allele count for inclusion. Default: 0.} +} +\value{ +3 x 1 numeric vector with the counts of BMVs, DMVs, and PTVs. +} +\description{ +Count the number of non-zero alleles bearing each variant annotation. +} diff --git a/tests/testthat/test-allelic_series.R b/tests/testthat/test-allelic_series.R index 4462b05..99c7fab 100644 --- a/tests/testthat/test-allelic_series.R +++ b/tests/testthat/test-allelic_series.R @@ -193,8 +193,8 @@ test_that("Overall check of omnibus test.", { include_orig_skato_all = FALSE, include_orig_skato_ptv = FALSE ) - expect_equal(length(p_omni), 8) - expect_equal(p_omni["p_omni"], 0.0, ignore_attr = TRUE, tolerance = 0.005) + expect_equal(length(p_omni), 3 + 8) + expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) # With SKAT-O all. p_omni <- COAST( @@ -204,8 +204,8 @@ test_that("Overall check of omnibus test.", { include_orig_skato_all = TRUE, include_orig_skato_ptv = FALSE ) - expect_equal(length(p_omni), 9) - expect_equal(p_omni["p_omni"], 0.0, ignore_attr = TRUE, tolerance = 0.005) + expect_equal(length(p_omni), 3 + 8 + 1) + expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) # With SKAT-O PTV. p_omni <- COAST( @@ -215,8 +215,8 @@ test_that("Overall check of omnibus test.", { include_orig_skato_all = FALSE, include_orig_skato_ptv = TRUE ) - expect_equal(length(p_omni), 9) - expect_equal(p_omni["p_omni"], 0.0, ignore_attr = TRUE, tolerance = 0.005) + expect_equal(length(p_omni), 3 + 8 + 1) + expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) # With both SKAT-O all and SKAT-O PTV. p_omni <- COAST( @@ -226,8 +226,8 @@ test_that("Overall check of omnibus test.", { include_orig_skato_all = TRUE, include_orig_skato_ptv = TRUE ) - expect_equal(length(p_omni), 10) - expect_equal(p_omni["p_omni"], 0.0, ignore_attr = TRUE, tolerance = 0.005) + expect_equal(length(p_omni), 3 + 8 + 2) + expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) }) diff --git a/tests/testthat/test-binary.R b/tests/testthat/test-binary.R index 4989aea..f749f6c 100644 --- a/tests/testthat/test-binary.R +++ b/tests/testthat/test-binary.R @@ -22,7 +22,7 @@ test_that("Overall check for binary phenotype.", { include_orig_skato_ptv = TRUE ) )) - expect_gt(p_omni["p_omni"], 0.05) + expect_gt(as.numeric(p_omni["p_omni"]), 0.05) # Alternative phenotype. invisible(capture.output( @@ -35,6 +35,6 @@ test_that("Overall check for binary phenotype.", { include_orig_skato_ptv = TRUE ) )) - expect_equal(p_omni["p_omni"], 0.0, ignore_attr = TRUE, tolerance = 0.005) + expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) }) \ No newline at end of file diff --git a/tests/testthat/test-counts.R b/tests/testthat/test-counts.R new file mode 100644 index 0000000..d797677 --- /dev/null +++ b/tests/testthat/test-counts.R @@ -0,0 +1,44 @@ +test_that("Check variant counter.", { + + # Data. + anno <- c(0, 1, 2) + geno <- rbind( + c(1, 0, 0), + c(0, 1, 0), + c(0, 2, 0), + c(0, 0, 1), + c(0, 0, 1), + c(0, 0, 2) + ) + + # No MAC filter: count alleles. + obs <- CountAlleles(anno, geno, count_carriers = FALSE, min_mac = 0) + exp <- c(n_bmv = 1, n_dmv = 3, n_ptv = 4) + expect_equal(obs, exp) + + # No MAC filter: count carriers. + obs <- CountAlleles(anno, geno, count_carriers = TRUE, min_mac = 0) + exp <- c(n_bmv = 1, n_dmv = 2, n_ptv = 3) + expect_equal(obs, exp) + + # MAC filter of 2: count alleles. + obs <- CountAlleles(anno, geno, count_carriers = FALSE, min_mac = 2) + exp <- c(n_bmv = 0, n_dmv = 3, n_ptv = 4) + expect_equal(obs, exp) + + # MAC filter of 2: count carriers. + obs <- CountAlleles(anno, geno, count_carriers = TRUE, min_mac = 2) + exp <- c(n_bmv = 0, n_dmv = 2, n_ptv = 3) + expect_equal(obs, exp) + + # MAC filter of 100: count alleles. + obs <- CountAlleles(anno, geno, count_carriers = FALSE, min_mac = 100) + exp <- c(n_bmv = 0, n_dmv = 0, n_ptv = 0) + expect_equal(obs, exp) + + # MAC filter of 100: count carriers. + obs <- CountAlleles(anno, geno, count_carriers = TRUE, min_mac = 100) + exp <- c(n_bmv = 0, n_dmv = 0, n_ptv = 0) + expect_equal(obs, exp) + +}) diff --git a/vignettes/coast.Rmd b/vignettes/coast.Rmd index 1fdff02..1c36499 100644 --- a/vignettes/coast.Rmd +++ b/vignettes/coast.Rmd @@ -1,7 +1,7 @@ --- title: "Coding-variant Allelic Series Test" output: rmarkdown::html_vignette -date: "2023-05-01" +date: "2024-04-03" vignette: > %\VignetteIndexEntry{my-vignette} %\VignetteEngine{knitr::rmarkdown} @@ -58,7 +58,7 @@ The example data generated by the preceding are available under `vignettes/vigne ## Running the alleic series test -The COding-variant Allelic Series Test (COAST) is run using the `COAST` function. By default, p-values for the component tests, as well as the overall omnibus test (`p_omni`), are returned. Inspection of the component p-values is useful for determining which model(s) drove an association. In the presence case, the association was most evident via the baseline count model (`p_count`). +The COding-variant Allelic Series Test (COAST) is run using the `COAST` function. By default, counts for the number of alleles of each variant class that contributed to the test are returned, along with p-values for the component tests, as well as the overall omnibus test (`p_omni`). Inspection of the component p-values may be useful for determining which model(s) detected evidence of an association. In the present case, the baseline count model (`p_count`) provided the greatest power. ```{r} results <- AllelicSeries::COAST( @@ -109,6 +109,17 @@ AllelicSeries::COAST( ) ``` +* `min_mac` is used to filter the variant set to those containing a minimum minor allele count (MAC). The following example filters to only those variants with a MAC of at least 2: +```{r, eval = FALSE} +AllelicSeries::COAST( + anno = anno, + geno = geno, + pheno = pheno, + covar = covar, + min_mac = 2 +) +``` + * `return_omni_only = TRUE` is used to return `p_omni` only when the component p-values are not of interest: ```{r, eval = FALSE}