Class for handling and visualizing results generated using rvat.

The rvatResult class is specifically designed to represent association results generated using RVAT (with assocTest or geneSetAssoc for example). It extends the BioConductor DFrame class, and allows for basic operations such as subsetting and merging as well as visualization (manhattan, qqplot, forest plots) and downstream analyses (e.g. ACAT). Different type of results have their own subclasses (rvbResult, singlevarResult, gsaResult) that inherit from rvatResult.

Constructing

rvbResult(object, header = TRUE): Intitialize an rvbResult object, including singlevar results as generated by the assocTest method. object can be either a data.frame/DataFrame or a filepath pointing to the results.
singlevarResult(object, header = TRUE): Intitialize an singlevarResult object, including rvb results as generated by the assocTest method. object can be either a data.frame/DataFrame or a filepath pointing to the results.
gsaResult(object, header = TRUE): Intitialize an singlevarResult object, including rvb results as generated by the geneSetAssoc or assocTest methods. object can be either a data.frame/DataFrame or a filepath pointing to the results
readResults(path, header = TRUE, type = NULL, sep = "\t"): Alternatively, results can be read using readResults where type can be set to rvbResult/singlevarResult/gsaResult. If type is not set, it will be inferred from the input.

Subsetting

Subsetting an rvatResult object is equivalent to subsetting a DataFrame/DFrame object and is fully described in S4Vectors::DataFrame. DataFrame objects behave very similar to the base data.frame, with the most notable exceptions that row names are optional and it can hold alternative vectors such as run-length encoded vectors (S4Vectors::Rle). Please see the example section for examples.

Merging

merge(x, y, by): Merge an rvatResult object with a data.frame or DataFrame object. The by parameter specifies which variables to join by. For example, by = c("a" = "b") will match x$a to y$b. Join by multiple variable by providing a vector of length>1 to by.

Displaying

In the following code snippets, x is an rvatResult object

show(x): By the default the number of rows and columns will be displayed + the first five rows.
summary(x, asList = FALSE): Shows a summary of the contents of the rvatResult. Optionally the output can be stored as a list by setting asList = TRUE.
topResult(x, n=10): Show the top N most significant results (based on the P-value).

Visualization

qqplot: Plot a qqplot based on the P-values stored in an rvatResult object. See qqplot for details
manhattan: Generate a manhattan plot for an rvatResult object. See manhattan for details.
rvatViewer: Interactive visualization, see rvatViewer for details.

Downsteam analyses

geneSetAssoc: An rvbResult can be used as input for geneSetAssoc to perform gene set analyses. See geneSetAssoc for details.
ACAT: P-values in an rvbResult can be combined using the ACAT method. See ACAT for details.

Writing

writeResult(x, ...): Write an rvatResult to disk, see writeResult for details.

Getters

getGdbId(x): Get the gdb ID
getGenomeBuild(x): Get gdb genome build

Examples


library(rvatData)
data(rvbresults)

# rvatResult inherits from the DataFrame classes
# standard methods also work on an rvatResult
head(rvbresults)
#> rvbResult with 6 rows and 28 columns
#>          unit cohort     varSetName   name pheno           covar geneticModel
#>   <character>  <Rle>          <Rle>  <Rle> <Rle>           <Rle>        <Rle>
#> 1        A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 2        A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 3        A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 4        A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 5        A1BG  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 6        A1BG  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#>   MAFweight               test      nvar caseCarriers ctrlCarriers
#>       <Rle>              <Rle> <numeric>    <numeric>    <numeric>
#> 1         1              firth        74           87          333
#> 2         1        skat_robust        74           87          333
#> 3         1           acatvSPA        74           87          333
#> 4         1 skat_burden_robust        74           87          333
#> 5         1              firth         4            2           13
#> 6         1        skat_robust         4            2           13
#>   meanCaseScore meanCtrlScore     caseN     ctrlN caseCallRate ctrlCallRate
#>       <numeric>     <numeric> <numeric> <numeric>    <numeric>    <numeric>
#> 1   0.017500262    0.01675191      5000     20000     0.993846     0.993780
#> 2   0.017500262    0.01675191      5000     20000     0.993846     0.993780
#> 3   0.017500262    0.01675191      5000     20000     0.993846     0.993780
#> 4   0.017500262    0.01675191      5000     20000     0.993846     0.993780
#> 5   0.000404543    0.00065511      5000     20000     0.995350     0.994775
#> 6   0.000404543    0.00065511      5000     20000     0.995350     0.994775
#>       effect  effectSE effectCIlower effectCIupper        OR         P
#>    <numeric> <numeric>     <numeric>     <numeric> <numeric> <numeric>
#> 1  0.0541596        NA     -0.189003     -1.249277  1.055653  0.656635
#> 2         NA        NA            NA            NA        NA  0.385580
#> 3         NA        NA            NA            NA        NA  0.446330
#> 4         NA        NA            NA            NA        NA  0.681784
#> 5 -0.3175556        NA     -1.945522     -0.133086  0.727926  0.632891
#> 6         NA        NA            NA            NA        NA  0.884041
#>       CHROM       POS     start       end
#>   <numeric> <numeric> <numeric> <numeric>
#> 1        19  58349335  58345178  58353492
#> 2        19  58349335  58345178  58353492
#> 3        19  58349335  58345178  58353492
#> 4        19  58349335  58345178  58353492
#> 5        19  58349335  58345178  58353492
#> 6        19  58349335  58345178  58353492
nrow(rvbresults)
#> [1] 139740
ncol(rvbresults)
#> [1] 28
dim(rvbresults)
#> [1] 139740     28
rvbresults_df <- as.data.frame(rvbresults)
head(rvbresults_df)
#>   unit cohort     varSetName   name pheno           covar geneticModel
#> 1 A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 2 A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 3 A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 4 A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 5 A1BG  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 6 A1BG  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#>   MAFweight               test nvar caseCarriers ctrlCarriers meanCaseScore
#> 1         1              firth   74           87          333  0.0175002617
#> 2         1        skat_robust   74           87          333  0.0175002617
#> 3         1           acatvSPA   74           87          333  0.0175002617
#> 4         1 skat_burden_robust   74           87          333  0.0175002617
#> 5         1              firth    4            2           13  0.0004045432
#> 6         1        skat_robust    4            2           13  0.0004045432
#>   meanCtrlScore caseN ctrlN caseCallRate ctrlCallRate      effect effectSE
#> 1  0.0167519094  5000 20000    0.9938459    0.9937804  0.05415961       NA
#> 2  0.0167519094  5000 20000    0.9938459    0.9937804          NA       NA
#> 3  0.0167519094  5000 20000    0.9938459    0.9937804          NA       NA
#> 4  0.0167519094  5000 20000    0.9938459    0.9937804          NA       NA
#> 5  0.0006551104  5000 20000    0.9953500    0.9947750 -0.31755562       NA
#> 6  0.0006551104  5000 20000    0.9953500    0.9947750          NA       NA
#>   effectCIlower effectCIupper        OR         P CHROM      POS    start
#> 1    -0.1890029    -1.2492774 1.0556531 0.6566347    19 58349335 58345178
#> 2            NA            NA        NA 0.3855803    19 58349335 58345178
#> 3            NA            NA        NA 0.4463299    19 58349335 58345178
#> 4            NA            NA        NA 0.6817840    19 58349335 58345178
#> 5    -1.9455216    -0.1330857 0.7279262 0.6328908    19 58349335 58345178
#> 6            NA            NA        NA 0.8840408    19 58349335 58345178
#>        end
#> 1 58353492
#> 2 58353492
#> 3 58353492
#> 4 58353492
#> 5 58353492
#> 6 58353492
rvbresults[1:3,]
#> rvbResult with 3 rows and 28 columns
#>          unit cohort     varSetName   name pheno           covar geneticModel
#>   <character>  <Rle>          <Rle>  <Rle> <Rle>           <Rle>        <Rle>
#> 1        A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 2        A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 3        A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#>   MAFweight        test      nvar caseCarriers ctrlCarriers meanCaseScore
#>       <Rle>       <Rle> <numeric>    <numeric>    <numeric>     <numeric>
#> 1         1       firth        74           87          333     0.0175003
#> 2         1 skat_robust        74           87          333     0.0175003
#> 3         1    acatvSPA        74           87          333     0.0175003
#>   meanCtrlScore     caseN     ctrlN caseCallRate ctrlCallRate    effect
#>       <numeric> <numeric> <numeric>    <numeric>    <numeric> <numeric>
#> 1     0.0167519      5000     20000     0.993846      0.99378 0.0541596
#> 2     0.0167519      5000     20000     0.993846      0.99378        NA
#> 3     0.0167519      5000     20000     0.993846      0.99378        NA
#>    effectSE effectCIlower effectCIupper        OR         P     CHROM       POS
#>   <numeric>     <numeric>     <numeric> <numeric> <numeric> <numeric> <numeric>
#> 1        NA     -0.189003      -1.24928   1.05565  0.656635        19  58349335
#> 2        NA            NA            NA        NA  0.385580        19  58349335
#> 3        NA            NA            NA        NA  0.446330        19  58349335
#>       start       end
#>   <numeric> <numeric>
#> 1  58345178  58353492
#> 2  58345178  58353492
#> 3  58345178  58353492
rvbresults[,1:10]
#> rvbResult with 139740 rows and 10 columns
#>               unit cohort     varSetName   name pheno           covar
#>        <character>  <Rle>          <Rle>  <Rle> <Rle>           <Rle>
#> 1             A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4
#> 2             A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4
#> 3             A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4
#> 4             A1BG  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4
#> 5             A1BG  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4
#> ...            ...    ...            ...    ...   ...             ...
#> 139736        ZZZ3  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4
#> 139737        ZZZ3  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4
#> 139738        ZZZ3  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4
#> 139739        ZZZ3  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4
#> 139740        ZZZ3  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4
#>        geneticModel MAFweight               test      nvar
#>               <Rle>     <Rle>              <Rle> <numeric>
#> 1           allelic         1              firth        74
#> 2           allelic         1        skat_robust        74
#> 3           allelic         1           acatvSPA        74
#> 4           allelic         1 skat_burden_robust        74
#> 5           allelic         1              firth         4
#> ...             ...       ...                ...       ...
#> 139736      allelic         1 skat_burden_robust       147
#> 139737      allelic         1              firth        10
#> 139738      allelic         1        skat_robust        10
#> 139739      allelic         1           acatvSPA        10
#> 139740      allelic         1 skat_burden_robust        10
rvbresults[["unit"]][1:5]
#> [1] "A1BG" "A1BG" "A1BG" "A1BG" "A1BG"

# write results
file <- tempfile()
writeResult(rvbresults, file = file)
rvbresults <- rvbResult(file)

# similarly, for single variant results
data(GTsmall)
sv <- assocTest(
  GTsmall,
  covar = c("sex", paste0("PC", 1:4)),
  pheno = "pheno",
  test = "scoreSPA",
  singlevar = TRUE,
  verbose = FALSE
)
svresultfile <- tempfile()
writeResult(sv, file = svresultfile)
sv <- singlevarResult(sv)
head(sv)
#> singlevarResult with 6 rows and 24 columns
#>        VAR_id cohort varSetName  name pheno               covar geneticModel
#>   <character>  <Rle>      <Rle> <Rle> <Rle>               <Rle>        <Rle>
#> 1        1268  pheno    unnamed  none pheno sex,PC1,PC2,PC3,PC4      allelic
#> 2        1270  pheno    unnamed  none pheno sex,PC1,PC2,PC3,PC4      allelic
#> 3        1271  pheno    unnamed  none pheno sex,PC1,PC2,PC3,PC4      allelic
#> 4        1273  pheno    unnamed  none pheno sex,PC1,PC2,PC3,PC4      allelic
#> 5        1275  pheno    unnamed  none pheno sex,PC1,PC2,PC3,PC4      allelic
#> 6        1276  pheno    unnamed  none pheno sex,PC1,PC2,PC3,PC4      allelic
#>       test   caseMAC   ctrlMAC     caseMAF   ctrlMAF     caseN     ctrlN
#>      <Rle> <numeric> <numeric>   <numeric> <numeric> <numeric> <numeric>
#> 1 scoreSPA         2         0 0.001050420         0       952      3801
#> 2 scoreSPA         0         0 0.000000000         0       992      3955
#> 3 scoreSPA         0         0 0.000000000         0       884      3551
#> 4 scoreSPA         1         0 0.000500501         0       999      3954
#> 5 scoreSPA         0         0 0.000000000         0      1002      3998
#> 6 scoreSPA         0         0 0.000000000         0       967      3895
#>   caseCallRate ctrlCallRate effectAllele otherAllele    effect  effectSE
#>      <numeric>    <numeric>        <Rle>       <Rle> <numeric> <numeric>
#> 1     0.950100     0.950725            T           C        NA        NA
#> 2     0.990020     0.989245            C           G        NA        NA
#> 3     0.882236     0.888194            G        GCAT        NA        NA
#> 4     0.997006     0.988994            G           A        NA        NA
#> 5     1.000000     1.000000            T           A        NA        NA
#> 6     0.965070     0.974237            C           G        NA        NA
#>   effectCIlower effectCIupper        OR         P
#>       <numeric>     <numeric> <numeric> <numeric>
#> 1            NA            NA        NA 0.0252463
#> 2            NA            NA        NA        NA
#> 3            NA            NA        NA        NA
#> 4            NA            NA        NA        NA
#> 5            NA            NA        NA        NA
#> 6            NA            NA        NA        NA

# merge 
merge <- merge(rvbresults[1:23], as.data.frame(rvbresults[,c(1, 23:28)]), by = "unit")
#> Rows have been added to the rvatResult object

# show summary
summary(rvbresults)
#> rvbResult object
#> -------------------
#> n units = 18385
#> N cases = 5000
#> N controls = 20000
#> names = maf001
#> varSets = ModerateImpact;HighImpact
#> tests = firth;skat_robust;acatvSPA;skat_burden_robust
#> covars = PC1,PC2,PC3,PC4
#> geneticModels = allelic
#> MAFweights = 1
#> pheno = pheno
#> cohorts = pheno

# show top results
topResult(rvbresults, n = 10)
#> rvbResult with 10 rows and 28 columns
#>           unit cohort     varSetName   name pheno           covar geneticModel
#>    <character>  <Rle>          <Rle>  <Rle> <Rle>           <Rle>        <Rle>
#> 1         SOD1  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 2         SOD1  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 3         SOD1  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 4         SOD1  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 5         NEK1  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 6         NEK1  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 7         NEK1  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 8         NEK1  pheno     HighImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 9         OPTN  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#> 10        OPTN  pheno ModerateImpact maf001 pheno PC1,PC2,PC3,PC4      allelic
#>    MAFweight               test      nvar caseCarriers ctrlCarriers
#>        <Rle>              <Rle> <numeric>    <numeric>    <numeric>
#> 1          1              firth        35           47           13
#> 2          1 skat_burden_robust        35           47           13
#> 3          1           acatvSPA        35           47           13
#> 4          1        skat_robust        35           47           13
#> 5          1 skat_burden_robust        33           37           29
#> 6          1              firth        33           37           29
#> 7          1           acatvSPA        33           37           29
#> 8          1        skat_robust        33           37           29
#> 9          1 skat_burden_robust        64           53          100
#> 10         1              firth        64           53          100
#>    meanCaseScore meanCtrlScore     caseN     ctrlN caseCallRate ctrlCallRate
#>        <numeric>     <numeric> <numeric> <numeric>    <numeric>    <numeric>
#> 1     0.00953410   0.000781703      5000     20000     0.968016     0.968674
#> 2     0.00953410   0.000781703      5000     20000     0.968016     0.968674
#> 3     0.00953410   0.000781703      5000     20000     0.968016     0.968674
#> 4     0.00953410   0.000781703      5000     20000     0.968016     0.968674
#> 5     0.00753214   0.001577446      5000     20000     0.961109     0.961045
#> 6     0.00753214   0.001577446      5000     20000     0.961109     0.961045
#> 7     0.00753214   0.001577446      5000     20000     0.961109     0.961045
#> 8     0.00753214   0.001577446      5000     20000     0.961109     0.961045
#> 9     0.01106932   0.005271168      5000     20000     0.958704     0.959369
#> 10    0.01106932   0.005271168      5000     20000     0.958704     0.959369
#>       effect  effectSE effectCIlower effectCIupper        OR           P
#>    <numeric> <numeric>     <numeric>     <numeric> <numeric>   <numeric>
#> 1   2.655538        NA       2.07918     1.1940081  14.23264 1.00000e-32
#> 2         NA        NA            NA            NA        NA 1.17824e-22
#> 3         NA        NA            NA            NA        NA 1.42552e-14
#> 4         NA        NA            NA            NA        NA 2.03265e-12
#> 5         NA        NA            NA            NA        NA 5.06659e-11
#> 6   1.633979        NA       1.15294     0.7534585   5.12423 8.30275e-11
#> 7         NA        NA            NA            NA        NA 1.75659e-06
#> 8         NA        NA            NA            NA        NA 2.03876e-06
#> 9         NA        NA            NA            NA        NA 7.06826e-06
#> 10  0.769394        NA       0.43423     0.0901051   2.15846 1.22572e-05
#>        CHROM       POS     start       end
#>    <integer> <numeric> <integer> <numeric>
#> 1         21  31664298  31659666  31668931
#> 2         21  31664298  31659666  31668931
#> 3         21  31664298  31659666  31668931
#> 4         21  31664298  31659666  31668931
#> 5          4 169491168 169369704 169612632
#> 6          4 169491168 169369704 169612632
#> 7          4 169491168 169369704 169612632
#> 8          4 169491168 169369704 169612632
#> 9         10  13118878  13099449  13138308
#> 10        10  13118878  13099449  13138308

# qqplot
man <- qqplot(rvbresults[rvbresults$varSetName == "ModerateImpact" & rvbresults$test == "firth",],
                 label = "unit")

# generate a manhattan plot
man <- manhattan(rvbresults[rvbresults$varSetName == "ModerateImpact" & rvbresults$test == "firth",],
                 label = "unit", 
                 contigs = "GRCh38")


# see ?ACAT and ?geneSetAssoc for downstream analyses on rvatResults