ACAT.RdCombine P-values in an rvbResult object using the ACAT method.
For details on the ACAT method see: (Liu et al., 2019).
rvatResult object
Variable to ACAT P-values across. For example, if aggregate = 'test', P-values across statistical tests will be combined using ACAT. A vector with multiple column names can be specified. For example, by specifying aggregate = c('test', 'varSetName'), P-values across statistical tests and varSets will be combined using ACAT. A list with multiple items can be specified, in which case step-wise ACAT is applied. For example, by specifying aggregate = list('test', 'varSetName'), first P-values across statistical tests are combined, and the resulting P-values are then combined across varSets.
Variables to group by.
For example, if group = c('unit', 'varSetName') and aggregate = 'test',
for each unit-varSetName combination, P-values across statistical tests are combined.
Defaults to c("unit", "cohort", "varSetName","name", "unit", "pheno", "covar", "geneticModel", "MAFweight", "test"),
i.e. all grouping variables in an rvbResult object.
The variable(s) specified in aggregate, are excluded from the grouping variables.
Should P-values that are exactly zero or one be fixed? Defaults to TRUE.
The method used for fixing the P-values can be specified using the fixpval_method parameter.
Method used to fix p-value (if fixpval = TRUE). Methods include:
'minmax' = P-values that are exactly 1 are replaced by the maximum value below 1 present in the results;
P-values that are exactly 0 are replaced by the minimum value above 0 in the results.
'manual' = Specify the replacement P-values using fixpval_maxP and fixpval_minP.
'Liu' = Method recommended by the author of the ACAT R package
(see: https://github.com/yaowuliu/ACAT),
P-values of 1 are replaced by 1-1/d, where d is the number of p-values combined by ACAT.
Since no recommendation for P-values of 0 is given, these are replaced with the value specified using fixpval_minP.
Replace P-values that are exactly 1 with this P-value if fixpval_method = 'manual'.
Replace P-values that are exactly 0 with this P-value if fixpval_method = 'manual'
or fixpval_method = 'Liu'.
Show warnings? Defaults to TRUE.
Liu, Y. et al. ACAT: A Fast and Powerful p Value Combination Method for Rare-Variant Analysis in Sequencing Studies. The American Journal of Human Genetics 104, 410–421 (2019).
library(rvatData)
# this will first ACAT P-values across statistical tests
# and then ACAT these P-values across varSets
data(rvbresults)
rvbresults <- rvbresults[1:1000,]
ACAT(
rvbresults,
aggregate = list("test", "varSetName")
)
#> Warning: 3 P-values are exactly 1, these are set to 0.999110017684605
#> rvbResult with 129 rows and 24 columns
#> unit cohort varSetName name pheno covar geneticModel
#> <character> <Rle> <Rle> <Rle> <Rle> <Rle> <Rle>
#> 1 A1BG pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 2 A1CF pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 3 A2M pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 4 A2ML1 pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 5 A3GALT2 pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> ... ... ... ... ... ... ... ...
#> 125 ACADM pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 126 ACADS pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 127 ACADSB pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 128 ACADVL pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 129 ACAN pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> MAFweight test nvar caseCarriers ctrlCarriers meanCaseScore
#> <Rle> <Rle> <numeric> <numeric> <numeric> <numeric>
#> 1 1 ACAT NA NA NA NA
#> 2 1 ACAT NA NA NA NA
#> 3 1 ACAT NA NA NA NA
#> 4 1 ACAT NA NA NA NA
#> 5 1 ACAT NA NA NA NA
#> ... ... ... ... ... ... ...
#> 125 1 ACAT NA NA NA NA
#> 126 1 ACAT NA NA NA NA
#> 127 1 ACAT NA NA NA NA
#> 128 1 ACAT NA NA NA NA
#> 129 1 ACAT 84 NA NA NA
#> meanCtrlScore caseN ctrlN caseCallRate ctrlCallRate effect
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1 NA 5000 20000 NA NA NA
#> 2 NA 5000 20000 NA NA NA
#> 3 NA 5000 20000 NA NA NA
#> 4 NA 5000 20000 NA NA NA
#> 5 NA 5000 20000 NA NA NA
#> ... ... ... ... ... ... ...
#> 125 NA 5000 20000 NA NA NA
#> 126 NA 5000 20000 NA NA NA
#> 127 NA 5000 20000 NA NA NA
#> 128 NA 5000 20000 NA NA NA
#> 129 NA 5000 20000 NA NA NA
#> effectSE effectCIlower effectCIupper OR P
#> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1 NA NA NA NA 0.638617
#> 2 NA NA NA NA 0.986629
#> 3 NA NA NA NA 0.435556
#> 4 NA NA NA NA 0.773626
#> 5 NA NA NA NA 0.771293
#> ... ... ... ... ... ...
#> 125 NA NA NA NA 0.803170
#> 126 NA NA NA NA 0.972401
#> 127 NA NA NA NA 0.702937
#> 128 NA NA NA NA 0.530125
#> 129 NA NA NA NA 0.384451
# alternatively, by providing a vector P-values across aggregates will be combined in one stage
ACAT(
rvbresults,
aggregate = c("test", "varSetName")
)
#> Warning: 3 P-values are exactly 1, these are set to 0.999110017684605
#> rvbResult with 129 rows and 24 columns
#> unit cohort varSetName name pheno covar geneticModel
#> <character> <Rle> <Rle> <Rle> <Rle> <Rle> <Rle>
#> 1 A1BG pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 2 A1CF pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 3 A2M pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 4 A2ML1 pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 5 A3GALT2 pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> ... ... ... ... ... ... ... ...
#> 125 ACADM pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 126 ACADS pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 127 ACADSB pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 128 ACADVL pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 129 ACAN pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> MAFweight test nvar caseCarriers ctrlCarriers meanCaseScore
#> <Rle> <Rle> <numeric> <numeric> <numeric> <numeric>
#> 1 1 ACAT NA NA NA NA
#> 2 1 ACAT NA NA NA NA
#> 3 1 ACAT NA NA NA NA
#> 4 1 ACAT NA NA NA NA
#> 5 1 ACAT NA NA NA NA
#> ... ... ... ... ... ... ...
#> 125 1 ACAT NA NA NA NA
#> 126 1 ACAT NA NA NA NA
#> 127 1 ACAT NA NA NA NA
#> 128 1 ACAT NA NA NA NA
#> 129 1 ACAT 84 67 312 0.0136886
#> meanCtrlScore caseN ctrlN caseCallRate ctrlCallRate effect
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1 NA 5000 20000 NA NA NA
#> 2 NA 5000 20000 NA NA NA
#> 3 NA 5000 20000 NA NA NA
#> 4 NA 5000 20000 NA NA NA
#> 5 NA 5000 20000 NA NA NA
#> ... ... ... ... ... ... ...
#> 125 NA 5000 20000 NA NA NA
#> 126 NA 5000 20000 NA NA NA
#> 127 NA 5000 20000 NA NA NA
#> 128 NA 5000 20000 NA NA NA
#> 129 0.0158356 5000 20000 0.993555 0.993661 NA
#> effectSE effectCIlower effectCIupper OR P
#> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1 NA NA NA NA 0.638617
#> 2 NA NA NA NA 0.986629
#> 3 NA NA NA NA 0.435556
#> 4 NA NA NA NA 0.773626
#> 5 NA NA NA NA 0.771293
#> ... ... ... ... ... ...
#> 125 NA NA NA NA 0.803170
#> 126 NA NA NA NA 0.972401
#> 127 NA NA NA NA 0.702937
#> 128 NA NA NA NA 0.530125
#> 129 NA NA NA NA 0.384451
# Input P-value shouldn't be exactly 0 or 1 (see: https://github.com/yaowuliu/ACAT).
# By default, P-values that are exactly 0 or 1 are reset (`fixpval = TRUE`) to the
# minimum P-value (>0) and maximum P-value (<1) in the results.
# Alternatives include:
# manual minmax values:
ACAT(
rvbresults,
aggregate = list("test", "varSetName"),
fixpval = TRUE,
fixpval_method = "manual",
fixpval_maxP = 0.9999,
fixpval_minP = 1e-32
)
#> Warning: 3 P-values are exactly 1, these are set to 0.9999
#> rvbResult with 129 rows and 24 columns
#> unit cohort varSetName name pheno covar geneticModel
#> <character> <Rle> <Rle> <Rle> <Rle> <Rle> <Rle>
#> 1 A1BG pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 2 A1CF pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 3 A2M pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 4 A2ML1 pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 5 A3GALT2 pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> ... ... ... ... ... ... ... ...
#> 125 ACADM pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 126 ACADS pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 127 ACADSB pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 128 ACADVL pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 129 ACAN pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> MAFweight test nvar caseCarriers ctrlCarriers meanCaseScore
#> <Rle> <Rle> <numeric> <numeric> <numeric> <numeric>
#> 1 1 ACAT NA NA NA NA
#> 2 1 ACAT NA NA NA NA
#> 3 1 ACAT NA NA NA NA
#> 4 1 ACAT NA NA NA NA
#> 5 1 ACAT NA NA NA NA
#> ... ... ... ... ... ... ...
#> 125 1 ACAT NA NA NA NA
#> 126 1 ACAT NA NA NA NA
#> 127 1 ACAT NA NA NA NA
#> 128 1 ACAT NA NA NA NA
#> 129 1 ACAT 84 NA NA NA
#> meanCtrlScore caseN ctrlN caseCallRate ctrlCallRate effect
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1 NA 5000 20000 NA NA NA
#> 2 NA 5000 20000 NA NA NA
#> 3 NA 5000 20000 NA NA NA
#> 4 NA 5000 20000 NA NA NA
#> 5 NA 5000 20000 NA NA NA
#> ... ... ... ... ... ... ...
#> 125 NA 5000 20000 NA NA NA
#> 126 NA 5000 20000 NA NA NA
#> 127 NA 5000 20000 NA NA NA
#> 128 NA 5000 20000 NA NA NA
#> 129 NA 5000 20000 NA NA NA
#> effectSE effectCIlower effectCIupper OR P
#> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1 NA NA NA NA 0.638617
#> 2 NA NA NA NA 0.986629
#> 3 NA NA NA NA 0.435556
#> 4 NA NA NA NA 0.773626
#> 5 NA NA NA NA 0.771293
#> ... ... ... ... ... ...
#> 125 NA NA NA NA 0.803170
#> 126 NA NA NA NA 0.972401
#> 127 NA NA NA NA 0.702937
#> 128 NA NA NA NA 0.530125
#> 129 NA NA NA NA 0.384451
# Liu method (see FAQ on: https://github.com/yaowuliu/ACAT)
ACAT(
rvbresults,
aggregate = list("test", "varSetName"),
fixpval = TRUE,
fixpval_method = "Liu"
)
#> rvbResult with 129 rows and 24 columns
#> unit cohort varSetName name pheno covar geneticModel
#> <character> <Rle> <Rle> <Rle> <Rle> <Rle> <Rle>
#> 1 A1BG pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 2 A1CF pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 3 A2M pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 4 A2ML1 pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 5 A3GALT2 pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> ... ... ... ... ... ... ... ...
#> 125 ACADM pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 126 ACADS pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 127 ACADSB pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 128 ACADVL pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> 129 ACAN pheno ACAT maf001 pheno PC1,PC2,PC3,PC4 allelic
#> MAFweight test nvar caseCarriers ctrlCarriers meanCaseScore
#> <Rle> <Rle> <numeric> <numeric> <numeric> <numeric>
#> 1 1 ACAT NA NA NA NA
#> 2 1 ACAT NA NA NA NA
#> 3 1 ACAT NA NA NA NA
#> 4 1 ACAT NA NA NA NA
#> 5 1 ACAT NA NA NA NA
#> ... ... ... ... ... ... ...
#> 125 1 ACAT NA NA NA NA
#> 126 1 ACAT NA NA NA NA
#> 127 1 ACAT NA NA NA NA
#> 128 1 ACAT NA NA NA NA
#> 129 1 ACAT 84 NA NA NA
#> meanCtrlScore caseN ctrlN caseCallRate ctrlCallRate effect
#> <numeric> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1 NA 5000 20000 NA NA NA
#> 2 NA 5000 20000 NA NA NA
#> 3 NA 5000 20000 NA NA NA
#> 4 NA 5000 20000 NA NA NA
#> 5 NA 5000 20000 NA NA NA
#> ... ... ... ... ... ... ...
#> 125 NA 5000 20000 NA NA NA
#> 126 NA 5000 20000 NA NA NA
#> 127 NA 5000 20000 NA NA NA
#> 128 NA 5000 20000 NA NA NA
#> 129 NA 5000 20000 NA NA NA
#> effectSE effectCIlower effectCIupper OR P
#> <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1 NA NA NA NA 0.638617
#> 2 NA NA NA NA 0.986629
#> 3 NA NA NA NA 0.435556
#> 4 NA NA NA NA 0.773626
#> 5 NA NA NA NA 0.771293
#> ... ... ... ... ... ...
#> 125 NA NA NA NA 0.803170
#> 126 NA NA NA NA 0.972401
#> 127 NA NA NA NA 0.702937
#> 128 NA NA NA NA 0.530125
#> 129 NA NA NA NA 0.384451