To run an allelic series test, there are 4 key inputs:
A numeric annotation vector, of the same length as the number of variants, coded as 0 for benign missense variants (BMVs), 1 for deleterious missense variants (DMVs), and 2 for protein truncating variants (PTVs).
A covariates matrix, with as many rows as subjects and including columns such as age and sex. If omitted, defaults to an intercept only.
A genotype matrix, with subjects as rows and variants as columns. The number of columns should correspond to the length of the annotation vector.
A numeric phenotype vector, either continuous or binary.
The example data used below were generated using the DGP
function provided with the package. The data set includes 100 subjects,
300 variants, and a continuous phenotype. The true effect sizes follow
an allelic series, with magnitudes proportional to
c(1, 2, 3) for BMVs, DMVs, and PTVs respectively.
set.seed(101)
n <- 100
data <- AllelicSeries::DGP(
n = n,
snps = 300,
beta = c(1, 2, 3) / sqrt(n),
)
# Annotations.
anno <- data$anno
head(anno)## [1] 0 0 0 2 1 0# Covariates.
covar <- data$covar
head(covar)## int age sex pc1 pc2 pc3
## [1,] 1 0.68477658 -1.1224972 -0.6807247 -0.4001798 -0.14686205
## [2,] 1 -0.39497721 0.8819621 0.4276220 0.4275115 -1.31046470
## [3,] 1 -0.23491732 0.8819621 2.6045440 1.3567071 -2.09212883
## [4,] 1 0.13831934 0.8819621 0.2811747 -0.3707890 0.07606076
## [5,] 1 -0.02007142 0.8819621 1.1725315 0.1239303 -1.46471403
## [6,] 1 1.76061280 0.8819621 -0.5597953 0.3463956 0.01463249# Genotypes.
geno <- data$geno
head(geno[,1:5])## [,1] [,2] [,3] [,4] [,5]
## [1,] 0 0 0 0 0
## [2,] 0 0 0 0 0
## [3,] 0 0 0 0 0
## [4,] 0 0 0 0 0
## [5,] 0 0 0 0 0
## [6,] 0 1 0 0 1# Phenotype.
pheno <- data$pheno
head(pheno)## [1] 2.02718287 -0.65067678 -0.62554778 0.04889603 -1.08599635 0.74506676The example data generated by the preceding are available under
vignettes/vignette_data.
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).
results <- AllelicSeries::COAST(
anno = anno,
geno = geno,
pheno = pheno,
covar = covar
)
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-28apply_int = TRUE applies the rank-based inverse normal
transformation from the 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.AllelicSeries::COAST(
anno = anno,
geno = geno,
pheno = pheno,
covar = covar,
apply_int = TRUE
)include_orig_skato_all = TRUE includes standard SKAT-O
applied to all variants as a component of the omnibus test, while
include_orig_skato_ptv = TRUE includes standard SKAT-O
applied to PTVs only. Including standard SKAT-O as a component of the
omnibus test can improve power to detect associations between the
phenotype and genes that may not be allelic series.AllelicSeries::COAST(
anno = anno,
geno = geno,
pheno = pheno,
covar = covar,
include_orig_skato_all = TRUE,
include_orig_skato_ptv = TRUE,
)is_pheno_binary = TRUE is required to indicate that the
supplied phenotype is binary, and should be analyzed using a logistic
regression model.AllelicSeries::COAST(
anno = anno,
geno = geno,
pheno = 1 * (pheno > 0),
covar = covar,
is_pheno_binary = TRUE
)return_omni_only = TRUE is used to return
p_omni only when the component p-values are not of
interest:AllelicSeries::COAST(
anno = anno,
geno = geno,
pheno = pheno,
covar = covar,
return_omni_only = TRUE
)weights specifies the relative phenotypic effects of
BMVs, DMVs, and PTVs. An increasing pattern such as the default setting
of weights = c(1, 2, 3) targets allelic series. Setting
weights = c(1, 1, 1) would target a genetic architecture
where all variants have equivalent expected magnitudes.AllelicSeries::COAST(
anno = anno,
geno = geno,
pheno = pheno,
covar = covar,
weights = c(1, 2, 3)
)