From 65eaeb563d881d10afa585837f06464f9c1249d7 Mon Sep 17 00:00:00 2001 From: zmccaw-insitro <91571107+zmccaw-insitro@users.noreply.github.com> Date: Thu, 25 Jul 2024 13:57:49 -0400 Subject: [PATCH] Update v0.6.0 (#20) * Update v0.5.0 * Fix vignette typo. * v0.6.0 * Review allelic_series_test.R Christoph's suggestion * Update test-counts.R --- DESCRIPTION | 3 +- NAMESPACE | 5 +- NEWS.md | 5 ++ R/RcppExports.R | 11 ++++ R/allelic_series_test.R | 53 +++++++++------- R/class.R | 56 +++++++++++++++++ R/utilities.R | 50 --------------- README.md | 91 ++++++++++++++++++++++------ man/COAST-class.Rd | 17 ++++++ man/COAST-method.Rd | 14 +++++ man/COAST.Rd | 3 - man/Counts.Rd | 21 +++++++ man/print.COAST.Rd | 16 +++++ src/RcppExports.cpp | 14 +++++ src/counting.cpp | 52 ++++++++++++++++ tests/testthat/test-allelic_series.R | 24 +++++--- tests/testthat/test-binary.R | 11 +++- tests/testthat/test-counts.R | 48 +++++++-------- tests/testthat/test-score_test.R | 8 ++- vignettes/coast.Rmd | 13 +++- 20 files changed, 379 insertions(+), 136 deletions(-) create mode 100644 R/class.R create mode 100644 man/COAST-class.Rd create mode 100644 man/COAST-method.Rd create mode 100644 man/Counts.Rd create mode 100644 man/print.COAST.Rd create mode 100644 src/counting.cpp diff --git a/DESCRIPTION b/DESCRIPTION index cfa805a..5dbf86d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: AllelicSeries Title: Allelic Series Test -Version: 0.0.5.0 +Version: 0.0.6.0 Authors@R: c(person(given = "Zachary", family = "McCaw", @@ -18,6 +18,7 @@ Description: Implementation of gene-level rare variant association tests targeti License: BSD_3_clause + file LICENSE Encoding: UTF-8 Imports: + methods, Rcpp, RNOmni, SKAT diff --git a/NAMESPACE b/NAMESPACE index 018120d..2d8592d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,12 +1,15 @@ # Generated by roxygen2: do not edit by hand +S3method(print,COAST) export(ASBT) export(ASKAT) export(Aggregator) export(COAST) export(Comparator) -export(CountAlleles) +export(Counts) export(DGP) export(OLS) +exportClasses(COAST) importFrom(Rcpp,sourceCpp) +importFrom(methods,show) useDynLib(AllelicSeries, .registration = TRUE) diff --git a/NEWS.md b/NEWS.md index e0cd170..9afb8d8 100644 --- a/NEWS.md +++ b/NEWS.md @@ -1,3 +1,8 @@ + +## Version 0.6.0 +* Replaced the `CountAlleles` function with a similar (but faster) `Counts` function that returns counts the total number of alleles, variants, and carriers by variant class. +* Updated the formatting of the results output by the main `COAST` function. + ## 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). diff --git a/R/RcppExports.R b/R/RcppExports.R index 076a1ed..5e33924 100644 --- a/R/RcppExports.R +++ b/R/RcppExports.R @@ -1,6 +1,17 @@ # Generated by using Rcpp::compileAttributes() -> do not edit by hand # Generator token: 10BE3573-1514-4C36-9D1C-5A225CD40393 +#' Count Variants and Carriers +#' +#' @param anno (snps x 1) annotation vector with values in c(0, 1, 2). +#' @param geno (n x snps) genotype matrix. +#' @param min_mac Minimum minor allele count for inclusion. Default: 0. +#' @return Data.frame of allele, variant, and carrier counts. +#' @export +Counts <- function(anno, geno, min_mac = 0L) { + .Call(`_AllelicSeries_Counts`, anno, geno, min_mac) +} + #' Ordinary Least Squares #' #' Fits the standard OLS model. diff --git a/R/allelic_series_test.R b/R/allelic_series_test.R index 95c0939..a69a11e 100644 --- a/R/allelic_series_test.R +++ b/R/allelic_series_test.R @@ -1,5 +1,5 @@ # Purpose: Allelic series test. -# Updated: 2024-04-03 +# Updated: 2024-07-24 # Default weights. DEFAULT_WEIGHTS <- c(1, 2, 3) @@ -346,7 +346,6 @@ ASKAT <- function( #' 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. @@ -376,7 +375,6 @@ COAST <- function( 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 @@ -442,12 +440,12 @@ COAST <- function( # Collect p-values. p_burden <- c( - p_count = p_count, - p_ind = p_ind, - p_max_count = p_max_count, - p_max_ind = p_max_ind, - p_sum_count = p_sum_count, - p_sum_ind = p_sum_ind + count = p_count, + ind = p_ind, + max_count = p_max_count, + max_ind = p_max_ind, + sum_count = p_sum_count, + sum_ind = p_sum_ind ) # SKAT tests. @@ -461,13 +459,13 @@ COAST <- function( return_null_model = TRUE, weights = weights ) - p_skat <- c(p_allelic_skat = skat_list$p) + p_skat <- c(allelic_skat = skat_list$p) null <- skat_list$null # Standard SKAT-O all. if (include_orig_skato_all) { orig_skato_all <- SKAT::SKAT(geno, null, method = "SKATO") - p_skat <- c(p_skat, p_orig_skat_all = orig_skato_all$p.value) + p_skat <- c(p_skat, orig_skat_all = orig_skato_all$p.value) } # Standard SKAT-O PTV. @@ -475,7 +473,7 @@ COAST <- function( ptv <- geno[, anno == 2, drop = FALSE] if (ncol(ptv) > 0) { orig_skato_ptv <- SKAT::SKAT(ptv, null, method = "SKATO") - p_skat <- c(p_skat, p_orig_skat_ptv = orig_skato_ptv$p.value) + p_skat <- c(p_skat, orig_skat_ptv = orig_skato_ptv$p.value) } } @@ -483,22 +481,33 @@ COAST <- function( n_burden <- length(p_burden) n_skat <- length(p_skat) - p_val = c(p_burden, p_skat) + pvals = c(p_burden, p_skat) omni_weights <- c(rep(1, n_burden), rep(n_burden / n_skat, n_skat)) - p_omni <- RNOmni::OmniP(p = p_val, w = omni_weights) + p_omni <- RNOmni::OmniP(p = pvals, w = omni_weights) - # Output. + # Only omnibus p-value requested. if (return_omni_only) { out <- c(p_omni = p_omni) - } 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) } + # Format p-values. + pvals <- c(pvals, omni = p_omni) + df_pvals <- data.frame( + test = names(pvals), + type = c(rep("burden", n_burden), rep("skat", n_skat), "omni"), + pval = as.numeric(pvals) + ) + + # Variant counts. + counts <- Counts(anno = anno, geno = geno, min_mac = min_mac) + + # Output. + out <- methods::new( + Class = "COAST", + Counts = counts, + Pvals = df_pvals + ) return(out) } diff --git a/R/class.R b/R/class.R new file mode 100644 index 0000000..dc0a6d8 --- /dev/null +++ b/R/class.R @@ -0,0 +1,56 @@ +# Purpose: Output class. +# Updated: 2024-07-24 + +#' Allelic Series Output Class +#' +#' @slot Counts Allele, variant, and carrier counts. +#' @slot Pvals Result p-values. +#' @name COAST-class +#' @rdname COAST-class +#' @exportClass COAST +setClass( + Class = "COAST", + representation = representation( + Counts = "data.frame", + Pvals = "data.frame" + ) +) + + +#' Print Method for COAST Object. +#' +#' Print method for objects of class \code{COAST}. +#' +#' @param x An object of class \code{COAST}. +#' @param ... Unused. +#' @export +print.COAST <- function(x, ...) { + + # Counts. + cat("Counts:\n") + show(x@Counts) + cat("\n\n") + + # P-values. + cat("P-values:\n") + pvals <- x@Pvals + pvals$pval <- formatC(pvals$pval, format = "e", digits = 2) + show(pvals) + cat("\n\n") +} + + +#' Show Method for COAST Object +#' +#' @param object An object of class \code{COAST}. +#' @rdname COAST-method +#' @importFrom methods show +setMethod( + f = "show", + signature = c(object = "COAST"), + definition = function(object) { + print.COAST(x = object) + } +) + + diff --git a/R/utilities.R b/R/utilities.R index dadb60f..00258cb 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -1,56 +1,6 @@ # Purpose: Utility functions. # 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 #' #' @param y (n x 1) continuous phenotype vector diff --git a/README.md b/README.md index 4026ff4..7f9bcae 100644 --- a/README.md +++ b/README.md @@ -1,3 +1,7 @@ +README +================ +2024-07-24 + # Allelic Series This package implements gene-level rare variant association tests @@ -93,19 +97,57 @@ genotype matrix, and the phenotype vector. show(results) ``` - ## 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 + ## Counts: + ## anno alleles variants carriers + ## 1 0 287 163 96 + ## 2 1 162 109 82 + ## 3 2 61 28 45 + ## + ## + ## P-values: + ## test type pval + ## 1 count burden 3.11e-26 + ## 2 ind burden 1.32e-09 + ## 3 max_count burden 3.08e-10 + ## 4 max_ind burden 5.37e-09 + ## 5 sum_count burden 1.66e-20 + ## 6 sum_ind burden 2.55e-11 + ## 7 allelic_skat skat 2.66e-07 + ## 8 omni omni 3.74e-25 + +By default, the output of `COAST` includes a data.frame of counts +showing the number of alleles, variants, and carriers in each class that +contributed to the test, and a data.frame of p-values, with the `omni` +test denoting the final, overall p-value. The counts data.frame is +accessed via: + +``` r +results@Counts +``` + + ## anno alleles variants carriers + ## 1 0 287 163 96 + ## 2 1 162 109 82 + ## 3 2 61 28 45 + +and the p-values data.frame via: + +``` r +results@Pvals +``` + + ## test type pval + ## 1 count burden 3.112702e-26 + ## 2 ind burden 1.322084e-09 + ## 3 max_count burden 3.076876e-10 + ## 4 max_ind burden 5.374363e-09 + ## 5 sum_count burden 1.661854e-20 + ## 6 sum_ind burden 2.554417e-11 + ## 7 allelic_skat skat 2.658137e-07 + ## 8 omni omni 3.735235e-25 -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`. +To return the omnibus p-value only, specify `return_omni_only = TRUE` +when calling `COAST`. ## Robust omnibus test @@ -133,12 +175,25 @@ results <- COAST( show(results) ``` - ## 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 + ## Counts: + ## anno alleles variants carriers + ## 1 0 287 163 96 + ## 2 1 162 109 82 + ## 3 2 61 28 45 + ## + ## + ## P-values: + ## test type pval + ## 1 count burden 3.11e-26 + ## 2 ind burden 1.32e-09 + ## 3 max_count burden 3.08e-10 + ## 4 max_ind burden 5.37e-09 + ## 5 sum_count burden 1.66e-20 + ## 6 sum_ind burden 2.55e-11 + ## 7 allelic_skat skat 2.66e-07 + ## 8 orig_skat_all skat 1.55e-05 + ## 9 orig_skat_ptv skat 6.63e-08 + ## 10 omni omni 3.74e-25 ## Loading genotypes diff --git a/man/COAST-class.Rd b/man/COAST-class.Rd new file mode 100644 index 0000000..05d7093 --- /dev/null +++ b/man/COAST-class.Rd @@ -0,0 +1,17 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/class.R +\docType{class} +\name{COAST-class} +\alias{COAST-class} +\title{Allelic Series Output Class} +\description{ +Allelic Series Output Class +} +\section{Slots}{ + +\describe{ +\item{\code{Counts}}{Allele, variant, and carrier counts.} + +\item{\code{Pvals}}{Result p-values.} +}} + diff --git a/man/COAST-method.Rd b/man/COAST-method.Rd new file mode 100644 index 0000000..8780bf3 --- /dev/null +++ b/man/COAST-method.Rd @@ -0,0 +1,14 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/class.R +\name{show,COAST-method} +\alias{show,COAST-method} +\title{Show Method for COAST Object} +\usage{ +\S4method{show}{COAST}(object) +} +\arguments{ +\item{object}{An object of class \code{COAST}.} +} +\description{ +Show Method for COAST Object +} diff --git a/man/COAST.Rd b/man/COAST.Rd index 37dc654..cb2ae4b 100644 --- a/man/COAST.Rd +++ b/man/COAST.Rd @@ -14,7 +14,6 @@ COAST( 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 @@ -42,8 +41,6 @@ to PTV variants only in the omnibus test? 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/Counts.Rd b/man/Counts.Rd new file mode 100644 index 0000000..3218821 --- /dev/null +++ b/man/Counts.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/RcppExports.R +\name{Counts} +\alias{Counts} +\title{Count Variants and Carriers} +\usage{ +Counts(anno, geno, min_mac = 0L) +} +\arguments{ +\item{anno}{(snps x 1) annotation vector with values in c(0, 1, 2).} + +\item{geno}{(n x snps) genotype matrix.} + +\item{min_mac}{Minimum minor allele count for inclusion. Default: 0.} +} +\value{ +Data.frame of allele, variant, and carrier counts. +} +\description{ +Count Variants and Carriers +} diff --git a/man/print.COAST.Rd b/man/print.COAST.Rd new file mode 100644 index 0000000..369884f --- /dev/null +++ b/man/print.COAST.Rd @@ -0,0 +1,16 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/class.R +\name{print.COAST} +\alias{print.COAST} +\title{Print Method for COAST Object.} +\usage{ +\method{print}{COAST}(x, ...) +} +\arguments{ +\item{x}{An object of class \code{COAST}.} + +\item{...}{Unused.} +} +\description{ +Print method for objects of class \code{COAST}. +} diff --git a/src/RcppExports.cpp b/src/RcppExports.cpp index 45af77f..f42cb9d 100644 --- a/src/RcppExports.cpp +++ b/src/RcppExports.cpp @@ -11,6 +11,19 @@ Rcpp::Rostream& Rcpp::Rcout = Rcpp::Rcpp_cout_get(); Rcpp::Rostream& Rcpp::Rcerr = Rcpp::Rcpp_cerr_get(); #endif +// Counts +SEXP Counts(arma::colvec anno, arma::mat geno, const int min_mac); +RcppExport SEXP _AllelicSeries_Counts(SEXP annoSEXP, SEXP genoSEXP, SEXP min_macSEXP) { +BEGIN_RCPP + Rcpp::RObject rcpp_result_gen; + Rcpp::RNGScope rcpp_rngScope_gen; + Rcpp::traits::input_parameter< arma::colvec >::type anno(annoSEXP); + Rcpp::traits::input_parameter< arma::mat >::type geno(genoSEXP); + Rcpp::traits::input_parameter< const int >::type min_mac(min_macSEXP); + rcpp_result_gen = Rcpp::wrap(Counts(anno, geno, min_mac)); + return rcpp_result_gen; +END_RCPP +} // OLS SEXP OLS(const arma::colvec y, const arma::mat X); RcppExport SEXP _AllelicSeries_OLS(SEXP ySEXP, SEXP XSEXP) { @@ -51,6 +64,7 @@ END_RCPP } static const R_CallMethodDef CallEntries[] = { + {"_AllelicSeries_Counts", (DL_FUNC) &_AllelicSeries_Counts, 3}, {"_AllelicSeries_OLS", (DL_FUNC) &_AllelicSeries_OLS, 2}, {"_AllelicSeries_ResidVar", (DL_FUNC) &_AllelicSeries_ResidVar, 2}, {"_AllelicSeries_Score", (DL_FUNC) &_AllelicSeries_Score, 4}, diff --git a/src/counting.cpp b/src/counting.cpp new file mode 100644 index 0000000..be0ff80 --- /dev/null +++ b/src/counting.cpp @@ -0,0 +1,52 @@ +// [[Rcpp::depends(RcppArmadillo)]] +#include + +// For debugging: Rcpp::Rcout << << std::endl; + +// ---------------------------------------------------------------------------- + +//' Count Variants and Carriers +//' +//' @param anno (snps x 1) annotation vector with values in c(0, 1, 2). +//' @param geno (n x snps) genotype matrix. +//' @param min_mac Minimum minor allele count for inclusion. Default: 0. +//' @return Data.frame of allele, variant, and carrier counts. +//' @export +// [[Rcpp::export]] +SEXP Counts( + arma::colvec anno, + arma::mat geno, + const int min_mac = 0 +){ + + // Apply MAC filter. + arma::rowvec macs = arma::sum(geno, 0); + geno = geno.cols(arma::find(macs.t() > min_mac)); + anno = anno.elem(arma::find(macs.t() > min_mac)); + + arma::colvec alleles = arma::zeros(3); + arma::colvec variants = arma::zeros(3); + arma::colvec carriers = arma::zeros(3); + +// Loop over annotations. + for(int j=0; j<3; j++) { + + arma::uvec key = arma::find(anno == j); + arma::mat geno_subset = geno.cols(key); + variants(j) = geno_subset.n_cols; + alleles(j) = arma::accu(geno_subset); + + // Total number of variants carried by each subject. + arma::colvec col_sum = arma::sum(geno_subset, 1); + carriers(j) = arma::accu(col_sum != 0); + + }; + + return Rcpp::DataFrame::create( + Rcpp::Named("anno")=arma::linspace(0, 2, 3), + Rcpp::Named("alleles")=alleles, + Rcpp::Named("variants")=variants, + Rcpp::Named("carriers")=carriers + ); +} + diff --git a/tests/testthat/test-allelic_series.R b/tests/testthat/test-allelic_series.R index 99c7fab..8acd871 100644 --- a/tests/testthat/test-allelic_series.R +++ b/tests/testthat/test-allelic_series.R @@ -193,8 +193,10 @@ test_that("Overall check of omnibus test.", { include_orig_skato_all = FALSE, include_orig_skato_ptv = FALSE ) - expect_equal(length(p_omni), 3 + 8) - expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) + pvals <- p_omni@Pvals + p_omni <- as.numeric(pvals$pval[pvals$test == "omni"]) + expect_equal(length(pvals$pval), 8) + expect_equal(p_omni, 0.0, tolerance = 0.005) # With SKAT-O all. p_omni <- COAST( @@ -204,8 +206,10 @@ test_that("Overall check of omnibus test.", { include_orig_skato_all = TRUE, include_orig_skato_ptv = FALSE ) - expect_equal(length(p_omni), 3 + 8 + 1) - expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) + pvals <- p_omni@Pvals + p_omni <- as.numeric(pvals$pval[pvals$test == "omni"]) + expect_equal(length(pvals$pval), 8 + 1) + expect_equal(p_omni, 0.0, tolerance = 0.005) # With SKAT-O PTV. p_omni <- COAST( @@ -215,8 +219,10 @@ test_that("Overall check of omnibus test.", { include_orig_skato_all = FALSE, include_orig_skato_ptv = TRUE ) - expect_equal(length(p_omni), 3 + 8 + 1) - expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) + pvals <- p_omni@Pvals + p_omni <- as.numeric(pvals$pval[pvals$test == "omni"]) + expect_equal(length(pvals$pval), 8 + 1) + expect_equal(p_omni, 0.0, tolerance = 0.005) # With both SKAT-O all and SKAT-O PTV. p_omni <- COAST( @@ -226,8 +232,10 @@ test_that("Overall check of omnibus test.", { include_orig_skato_all = TRUE, include_orig_skato_ptv = TRUE ) - expect_equal(length(p_omni), 3 + 8 + 2) - expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) + pvals <- p_omni@Pvals + p_omni <- as.numeric(pvals$pval[pvals$test == "omni"]) + expect_equal(length(pvals$pval), 8 + 2) + expect_equal(p_omni, 0.0, tolerance = 0.005) }) diff --git a/tests/testthat/test-binary.R b/tests/testthat/test-binary.R index f749f6c..8a6d4ff 100644 --- a/tests/testthat/test-binary.R +++ b/tests/testthat/test-binary.R @@ -22,7 +22,9 @@ test_that("Overall check for binary phenotype.", { include_orig_skato_ptv = TRUE ) )) - expect_gt(as.numeric(p_omni["p_omni"]), 0.05) + pvals <- p_omni@Pvals + p_omni <- as.numeric(pvals$pval[pvals$test == "omni"]) + expect_gt(p_omni, 0.05) # Alternative phenotype. invisible(capture.output( @@ -35,6 +37,9 @@ test_that("Overall check for binary phenotype.", { include_orig_skato_ptv = TRUE ) )) - expect_equal(as.numeric(p_omni["p_omni"]), 0.0, tolerance = 0.005) + pvals <- p_omni@Pvals + p_omni <- as.numeric(pvals$pval[pvals$test == "omni"]) + expect_equal(p_omni, 0.0, tolerance = 0.005) -}) \ No newline at end of file +}) + \ No newline at end of file diff --git a/tests/testthat/test-counts.R b/tests/testthat/test-counts.R index d797677..0d792b9 100644 --- a/tests/testthat/test-counts.R +++ b/tests/testthat/test-counts.R @@ -11,34 +11,30 @@ test_that("Check variant counter.", { 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) + # Counts with minimum MAC = 0. + counts <- Counts(anno = anno, geno = geno, min_mac = 0) + expect_equal(counts$alleles, c(1, 1 + 2, 1 + 1 + 2)) + expect_equal(counts$variants, c(1, 1, 1)) + expect_equal(counts$carriers, c(1, 2, 3)) - # 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) + # Counts with minimum MAC = 2. + counts <- Counts(anno = anno, geno = geno, min_mac = 2) + expect_equal(counts$alleles, c(0, 1 + 2, 1 + 1 + 2)) + expect_equal(counts$variants, c(0, 1, 1)) + expect_equal(counts$carriers, c(0, 2, 3)) - # 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) + # Counts with minimum MAC = 100. + counts <- Counts(anno = anno, geno = geno, min_mac = 100) + expect_equal(counts$alleles, c(0, 0, 0)) + expect_equal(counts$variants, c(0, 0, 0)) + expect_equal(counts$carriers, c(0, 0, 0)) - # 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) + # Check case of multiple variants per annotation category. + anno <- c(0, 1, 1, 2, 2, 2) + geno <- diag(6) # Identity matrix. + counts <- Counts(anno = anno, geno = geno, min_mac = 0) + expect_equal(counts$alleles, c(1, 2, 3)) + expect_equal(counts$variants, c(1, 2, 3)) + expect_equal(counts$carriers, c(1, 2, 3)) }) diff --git a/tests/testthat/test-score_test.R b/tests/testthat/test-score_test.R index 8bee7a5..1320891 100644 --- a/tests/testthat/test-score_test.R +++ b/tests/testthat/test-score_test.R @@ -52,7 +52,9 @@ test_that("Overall check of score test.", { include_orig_skato_ptv = FALSE, score_test = TRUE ) - expect_gt(p_omni_null["p_omni"], 0.05) + pvals <- p_omni_null@Pvals + p_omni <- as.numeric(pvals$pval[pvals$test == "omni"]) + expect_gt(p_omni, 0.05) # Associated phenotype. @@ -65,6 +67,8 @@ test_that("Overall check of score test.", { include_orig_skato_ptv = FALSE, score_test = TRUE ) - expect_lt(p_omni_alt["p_omni"], 0.05) + pvals <- p_omni_alt@Pvals + p_omni <- as.numeric(pvals$pval[pvals$test == "omni"]) + expect_lt(p_omni, 0.05) }) \ No newline at end of file diff --git a/vignettes/coast.Rmd b/vignettes/coast.Rmd index 1c36499..10b2070 100644 --- a/vignettes/coast.Rmd +++ b/vignettes/coast.Rmd @@ -1,7 +1,7 @@ --- title: "Coding-variant Allelic Series Test" output: rmarkdown::html_vignette -date: "2024-04-03" +date: "2024-07-24" 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, 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. +The COding-variant Allelic Series Test (COAST) is run using the `COAST` function. By default, the output of `COAST` includes a data.frame of counts showing the number of alleles, variants, and carriers in each class that contributed to the test, and a data.frame of p-values, with the `omni` test denoting the final, overall p-value. 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 provided the greatest power. ```{r} results <- AllelicSeries::COAST( @@ -70,6 +70,15 @@ results <- AllelicSeries::COAST( show(results) ``` +The counts data.frame is accessed via: +```{r} +results@Counts +``` +and the p-values data.frame via: +```{r} +results@Pvals +``` + ### Test options * `apply_int = TRUE` applies the rank-based inverse normal transformation from the [RNOmni](https://CRAN.R-project.org/package=RNOmni) package. This transformation is expected to improve power for phenotypes that have a skewed or kurtotic (e.g. long-tailed) distribution. It is applied by default in the case of continuous phenotype, and is ignored in the case of a binary phenotype.