assocTestAggregate {GENESIS} | R Documentation

## Aggregate Association Testing

### Description

`assocTestAggregate` performs aggregate association tests using the null model fit with `fitNullModel`.

### Usage

    ## S4 method for signature 'SeqVarIterator'
    assocTestAggregate(gdsobj, null.model, AF.max=1,
                       weight.beta=c(1,1), weight.user=NULL,
                       test=c("Burden", "SKAT", "SMMAT"),
                       burden.test=c("Score", "Wald"), rho=0,
                       pval.method=c("davies", "kuonen", "liu"),
                       sparse=TRUE, verbose=TRUE)
    ## S4 method for signature 'GenotypeIterator'
    assocTestAggregate(gdsobj, null.model, AF.max=1,
                       weight.beta=c(1,1), weight.user=NULL,
                       test=c("Burden", "SKAT", "SMMAT"),
                       burden.test=c("Score", "Wald"), rho=0,
                       pval.method=c("davies", "kuonen", "liu"),
                       verbose=TRUE)

### Arguments

`gdsobj` | An object of class `SeqVarIterator` from the package SeqVarTools containing the genotype data for the variants and samples to be used for the analysis.
`null.model` | A null model object returned by `fitNullModel`.
`AF.max` | A numeric value specifying the upper bound on the alternate allele frequency for variants to be included in the analysis.
`weight.beta` | A numeric vector of length two specifying the two parameters of the Beta distribution used to determine variant weights; weights are given by `dbeta(MAF, a, b)`, where MAF is the minor allele frequency, and a and b are the two parameters specified here. `weight.beta = c(1,25)` gives the Wu weights; `weight.beta = c(0.5, 0.5)` is proportional to the Madsen-Browning weights; and `weight.beta = c(1,1)` gives a weight of 1 to all variants. This input is ignored when `weight.user` is not `NULL`.
`weight.user` | A character string specifying the name of a variable to be used as variant weights. This variable can be in either 1) the variantData slot of `gdsobj` or 2) the `mcols` of the `GRanges` or `GRangesList` object used to create `gdsobj` (when `gdsobj` is a `link{SeqVarRangeIterator}` or `link{SeqVarListIterator}`). When left `NULL` (the default), the weights specified by `weight.beta` will be used.
`test` | A character string specifying the type of test to be performed. The possibilities are `"Burden"` (default), `"SKAT"`, or `"SMMAT"`. When this is set to "SKAT" and the parameter `rho` has multiple values, a SKAT-O test is performed.
`burden.test` | A character string specifying the type of Burden test to perform when `test = "Burden"`. The possibilities are `"Score"` and `"Wald"`. `"Score"` can be used for any `null.model`. `"Wald"` can not be used when the `null.model` is from a mixed model with a binary outcome variable.

`rho` | A numeric value (or vector of numeric values) in `[0,1]` specifying the rho parameter for SKAT. When `rho = 0`, a standard SKAT test is performed. When `rho = 1`, a score burden test is performed. When `rho` is a vector of values, SKAT-O is performed using each of those values as the search space for the optimal `rho`.
`pval.method` | A character string specifying which method to use to calculate SKAT p-values. `"davies"` (the default) uses numerical integration; `"kuonen"` uses a saddlepoint method; and `"liu"` uses a moment matching approximation. If the davies method generates an error, kuonen is tried, and then liu as a last resort.
`sparse` | Logical indicator of whether to read genotypes as sparse Matrix objects; the default is `TRUE`. Set this to `FALSE` if the alternate allele dosage of the genotypes in the test are not expected to be mostly 0.
`verbose` | Logical indicator of whether updates from the function should be printed to the console; the default is `TRUE`.

### Details

The type of aggregate unit tested depends on the class of iterator used for `gdsobj`. Options include sliding windows, specific ranges of variants or selection of individual variants (ranges with width 1). See `SeqVarIterator` for more details.

The effect size estimate is for each copy of the alternate allele. For multiallelic variants, each alternate allele is tested separately.

Somewhat similarly to SKAT-O, the variant Set Mixed Model Association Test (SMMAT, Chen et al., manuscript in preparation) combines the burden test p-value with an adjusted SKAT (which is asymptotically independent of the burden test) p-value using a chi-square distribution with 4df from Fisher's method.

### Value

A list with the following items:

`results` | A data.frame containing the results from the main analysis. Each row is a separate aggregate test:

If `gdsobj` is a `SeqVarWindowIterator`:

`chr` | The chromosome value
`start` | The start position of the window
`end` | The end position of the window

Always:

`n.site` | The number of variant sites included in the test.
`n.alt` | The number of alternate alleles included in the test.
`n.sample.alt` | The number of samples with an observed alternate allele at any variant in the aggregate set.

If `test` is `"Burden"`:

If `burden.test` is "Score":

`Score` | The value of the score function
`Score.SE` | The estimated standard error of the Score
`Score.Stat` | The score Z test statistic
`Score.pval` | The score p-value

If `burden.test` is `"Wald"`:

`Est` | The effect size estimate for a one unit increase in the burden value
`Est.SE` | The estimated standard error of the effect size estimate
`Wald.Stat` | The Wald Z test statistic
`Wald.pval` | The Wald p-value

If `test` is `"SKAT"`:

`Q_rho` | The SKAT test statistic for the value of rho specified. There will be as many of these variables as there are rho values chosen.
`pval_rho` | The SKAT p-value for the value of rho specified. There will be as many of these variables as there are rho values chosen.
`err_rho` | Takes value 1 if there was an error in calculating the p-value for the value of rho specified when using the "kuonen" or "davies" methods; 0 otherwise. When there is an error, the p-value returned is from the "liu" method. There will be as many of these variables as there are rho values chosen.

When `length(rho) > 1` and SKAT-O is performed:

`min.pval` | The minimum p-value among the p-values calculated for each choice of rho.
`opt.rho` | The optimal rho value; i.e. the rho value that gave the minimum p-value.
`pval_SKATO` | The SKAT-O p-value after adjustment for searching across multiple rho values.

If `test` is `"SMMAT"`:

`pval_burden` | The burden test p-value
`pval_SMMAT` | The SMMAT p-value
`err` | Takes value 1 if there was an error calculating the SMMAT p-value; 0 otherwise. If `err=1`, `pval_SMMAT` is set to `pval_burden`.
`variantInfo` | A list with as many elements as aggregate tests performed. Each element of the list is a data.frame providing information on the variants used in the aggregate test with results presented in the corresponding row of `results`. Each of these data.frames has the following information:
`variant.id` | The variant ID
`chr` | The chromosome value
`pos` | The base pair position
`n.obs` | The number of samples with non-missing genotypes
`freq` | The estimated alternate allele frequency
`weight` | The weight assigned to the variant in the analysis.

### Author(s)

Matthew P. Conomos, Stephanie M. Gogarten, Tamar Sofer, Ken Rice, Chaoyu Yu, Han Chen

### References

Leal, S.M. & Li, B. (2008). Methods for Detecting Associations with Rare Variants for Common Diseases: Application to Analysis of Sequence Data. American Journal of Human Genetics, 83(3): 311-321.

Browning, S.R. & Madsen, B.E. (2009). A Groupwise Association Test for Rare Mutations Using a Weighted Sum Statistic. PLoS Genetics, 5(2): e1000384.

Wu, M.C, Lee, S., Cai, T., Li, Y., Boehnke, M., & Lin, X. (2011). Rare-Variant Association Testing for Sequencing Data with the Sequence Kernel Association Test. American Journal of Human Genetics, 89(1): 82-93.

Lee, S. et al. (2012). Optimal Unified Approach for Rare-Variant Association Testing with Application to Small-Sample Case-Control Whole-Exome Sequencing Studies. American Journal of Human Genetics, 91(2): 224-237.

### Examples

    library(SeqVarTools)
    library(Biobase)
    library(GenomicRanges)

    # open a sequencing GDS file
    gdsfile <- seqExampleFileName("gds")
    gds <- seqOpen(gdsfile)

    # simulate some phenotype data
    data(pedigree)
    pedigree <- pedigree[match(seqGetData(gds, "sample.id"), pedigree$sample.id),]
    pedigree$outcome <- rnorm(nrow(pedigree))

    # construct a SeqVarData object
    seqData <- SeqVarData(gds, sampleData=AnnotatedDataFrame(pedigree))

    # fit the null model
    nullmod <- fitNullModel(seqData, outcome="outcome", covars="sex")

    # burden test - Range Iterator
    gr <- GRanges(seqnames=rep(1,3), ranges=IRanges(start=c(1e6, 2e6, 3e6), width=1e6))
    iterator <- SeqVarRangeIterator(seqData, variantRanges=gr)
    assoc <- assocTestAggregate(iterator, nullmod, test="Burden")
    assoc$results
    lapply(assoc$variantInfo, head)

    # SKAT test - Window Iterator
    seqSetFilterChrom(seqData, include="22")
    iterator <- SeqVarWindowIterator(seqData)
    assoc <- assocTestAggregate(iterator, nullmod, test="SKAT")
    head(assoc$results)
    head(assoc$variantInfo)

    # SKAT-O test - List Iterator
    seqResetFilter(iterator)
    gr <- GRangesList(
      GRanges(seqnames=rep(22,2), ranges=IRanges(start=c(16e6, 17e6), width=1e6)),
      GRanges(seqnames=rep(22,2), ranges=IRanges(start=c(18e6, 20e6), width=1e6)))
    iterator <- SeqVarListIterator(seqData, variantRanges=gr)
    assoc <- assocTestAggregate(iterator, nullmod, test="SKAT", rho=seq(0, 1, 0.25))
    assoc$results
    assoc$variantInfo

    # user-specified weights - option 1
    seqResetFilter(iterator)
    variant.id <- seqGetData(gds, "variant.id")
    weights <- data.frame(variant.id, weight=runif(length(variant.id)))
    variantData(seqData) <- AnnotatedDataFrame(weights)
    iterator <- SeqVarListIterator(seqData, variantRanges=gr)
    assoc <- assocTestAggregate(iterator, nullmod, test="Burden", weight.user="weight")
    assoc$results
    assoc$variantInfo

    # user-specified weights - option 2
    seqResetFilter(iterator)
    variantData(seqData)$weight <- NULL
    gr <- GRangesList(
      GRanges(seqnames=rep(22,2), ranges=IRanges(start=c(16e6, 17e6), width=1e6), weight=runif(2)),
      GRanges(seqnames=rep(22,2), ranges=IRanges(start=c(18e6, 20e6), width=1e6), weight=runif(2)))
    iterator <- SeqVarListIterator(seqData, variantRanges=gr)
    assoc <- assocTestAggregate(iterator, nullmod, test="Burden", weight.user="weight")
    assoc$results
    assoc$variantInfo

    seqClose(seqData)

[Package GENESIS version 2.12.4 Index]