Installation

QTLExperiment and multistateQTL can be installed from GitHub:

if (!require("BiocManager", quietly=TRUE))
    install.packages("BiocManager")

BiocManager::install(c("QTLExperiment", "multistateQTL"))

They are also available on GitHub:

devtools::install_git("https://github.com/dunstone-a/QTLExperiment", build_vignettes=TRUE)
devtools::install_git(https://github.com/dunstone-a/multistateQTL", build_vignettes=TRUE)

library(QTLExperiment)
library(multistateQTL)

Estimate parameters from GTEx

Provided with real QTL summary statistics as either a QTLE object or a named list with betas and error (and optionally pval or lfsr), key parameters are estimated that are used to simulate realistic multi-state QTL data. We demonstrate the parameter estimation function on publicly available summary statistics from GTEx (v8). Note that this data only contains tests that were called as significant by GTEx and vroom only loads the first chunk of results as it does not read .gz compressed objects well. This truncated dataset is used in the vignette for convenience, however, to estimate the default parameters in qtleParams() we downloaded QTL summary statistics for all associations tested for the 10 GTEx tissues with the largest sample sizes from Google Cloud). To speed up calculations we filtered this to include only associations on chromosome 1 and considered significant tests with pval < 0.05 and null tests with a pval > 0.1.

See QTLExperiment for more info on the sumstats2qtle function and for other approaches for reading in QTL summary statistics.

The parameters estimated here include:

• Significant betas shape and rate: Define the gamma distribution used to sample mean effect sizes for each QTL that is significant in at least one state. These simulated mean effect sizes are used as the mean parameter in rnorm to sample an effect size for each QTL for each state. The variance parameter for rnorm is user defined (default = 0.1).
  
• Significant coefficient of variation (cv) shape and rate: Define the gamma distribution used to sample the cv for each QTL in each state where that QTL is significant. The cv is multiplied by the simulated significant effect size for that test/state to get the simulated standard error values.
• Null beta shape and rate: Define the gamma distribution used to sample the effect sizes for each QTL in each state where the effect is not significant.
• Null beta cv shape and rate: Define the gamma distribution used to sample the cv for each QTL in each state where that QTL is not significant. This cv is then multiplied by the simulated null effect size for that test/state to get the simulated stand error values.

input_path <- system.file("extdata", package="multistateQTL")
state <- c("lung", "thyroid", "spleen", "blood")

input <- data.frame(
    state=state, 
    path=paste0(input_path, "/GTEx_tx_", state, ".tsv"))

gtex <- sumstats2qtle(
    input, 
    feature_id="molecular_trait_id",
    variant_id="rsid", 
    betas="beta", 
    errors="se",
    pvalues="pvalue", 
    verbose=TRUE)
gtex

## class: QTLExperiment 
## dim: 1163 4 
## metadata(0):
## assays(3): betas errors pvalues
## rownames(1163): ENST00000428771|rs554008981 ENST00000477976|rs554008981
##   ... ENST00000445118|rs368254419 ENST00000483767|rs368254419
## rowData names(2): feature_id variant_id
## colnames(4): lung thyroid spleen blood
## colData names(1): state_id

head(betas(gtex))

##                                   lung    thyroid   spleen     blood
## ENST00000428771|rs554008981 -0.1733690         NA 0.134913        NA
## ENST00000477976|rs554008981  0.1616170  0.3173110       NA        NA
## ENST00000483767|rs554008981 -0.4161480 -0.0483018       NA -0.204647
## ENST00000623070|rs554008981 -0.1137930         NA       NA        NA
## ENST00000669922|rs554008981 -0.1921680 -0.1067540 0.724622 -0.117424
## ENST00000428771|rs201055865 -0.0630909         NA       NA        NA

Estimating parameters:

params <- qtleEstimate(gtex, threshSig=0.05, threshNull=0.5)
params

## $cv.sig.shape
## [1] 7.070354
## 
## $cv.sig.rate
## [1] 7.717703
## 
## $cv.null.shape
## [1] 1.903592
## 
## $cv.null.rate
## [1] 0.8979238
## 
## $betas.sig.shape
## [1] 3.470245
## 
## $betas.sig.rate
## [1] 16.27222
## 
## $betas.null.shape
## [1] 1.360183
## 
## $betas.null.rate
## [1] 11.03394

Looking at the distributions defined by these estimated parameters, the simulated effect sizes for significant QTL will tend to be larger, while the simulated coefficient of variation values will be smaller than for the non-significant QTL.

plotSimulationParams(params=params)

Gamma distributions defined by the parameters estimated by qtleEstimate.

The default parameters available through qtleParams() were estimated from the GTEx v8 tissue-level eQTL summary statistics from chromosome 1 using the 10 tissues with the largest sample sizes. From these data, significant QTL parameters were estimated from tests in the lowest p-value quantile, while null parameters were estimated from tests in the highest p-value quantile. Data for tests on chromosome 1 were included in all four tissues (n=32613).

Simulate multi-state QTL data

The simulation tool allows for the simulation of four types of associations: (1) Global, where the simulated effect size is approximately equal across all states; (2) Unique, where the association is only significant in one state; (3) Multi-state, where the association is significant in a subset of states (i.e., state-groups), and (4) Null, where the association has no significant effects in any state. First each test is randomly assigned as one of the above types according to the proportions specified by the user. For multi-state QTL, each state is assigned to a state-group, either randomly or according to user defined groups, then each multi-state QTL is assigned randomly to one of the state-groups. For unique QTL, the QTL is randomly assigned to a single state.

Simulated mean effect sizes for all non-null QTL are sampled from gamma(beta.sig.shape, beta.sig.rate) and are randomly assigned a positive or negative effect direction. Then for each state where that QTL is significant, an effect size is sampled from N(mean effect size, σ), where σ is user defined (default=0.1). Effect sizes for null QTL are sampled from gamma(beta.null.shape, beta.null.rate) and are randomly assigned a positive or negative effect direction. Standard errors for each QTL for each state are simulated by sampling from gamma(cv.sig.shape, cv.sig.rate) or gamma(cv.null.shape, cv.null.rate) for significant and null QTL, respectively, and multiplying the sampled cv by the absolute value of the simulated beta for that QTL in that state.

Here is an example of a simple simulation with half of the simulated QTL tests having globally significant effects. This example uses the default parameters.

sim <- qtleSimulate(nTests=1000, nStates=6, global=0.5)
sim

## class: QTLExperiment 
## dim: 1000 6 
## metadata(0):
## assays(3): betas errors lfsrs
## rownames(1000): F0867|v86973 F0107|v67504 ... F0874|v42722 F0190|v78203
## rowData names(11): feature_id variant_id ... S5 S6
## colnames(6): S1 S2 ... S5 S6
## colData names(1): state_id

head(rowData(sim))

## DataFrame with 6 rows and 11 columns
##               feature_id  variant_id         QTL           id mean_beta
##              <character> <character> <character>  <character> <numeric>
## F0867|v86973       F0867      v86973      global F0867|v86973  0.274222
## F0107|v67504       F0107      v67504      global F0107|v67504  0.740930
## F0239|v64205       F0239      v64205      global F0239|v64205  0.211987
## F0531|v90951       F0531      v90951      global F0531|v90951 -0.650529
## F0950|v93596       F0950      v93596        null F0950|v93596  0.000000
## F0851|v25798       F0851      v25798        null F0851|v25798  0.000000
##                     S1        S2        S3        S4        S5        S6
##              <logical> <logical> <logical> <logical> <logical> <logical>
## F0867|v86973      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE
## F0107|v67504      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE
## F0239|v64205      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE
## F0531|v90951      TRUE      TRUE      TRUE      TRUE      TRUE      TRUE
## F0950|v93596     FALSE     FALSE     FALSE     FALSE     FALSE     FALSE
## F0851|v25798     FALSE     FALSE     FALSE     FALSE     FALSE     FALSE

We can also generate more complex simulations, for example this simulation has 20% global, 40% multi-state, 20% unique, and 20% null QTL effects, where multi-state effects are assigned to one of two state-groups.

sim <- qtleSimulate(
    nStates=10, nFeatures=100, nTests=1000,
    global=0.2, multi=0.4, unique=0.2, k=2)

Here is a snapshot of the simulation key for QTL simulated as unique to a single state:

head(rowData(subset(sim, QTL == "unique")))

## DataFrame with 6 rows and 16 columns
##              feature_id  variant_id         QTL          id mean_beta       S01
##             <character> <character> <character> <character> <numeric> <logical>
## F048|v14842        F048      v14842      unique F048|v14842 -0.159967     FALSE
## F087|v59649        F087      v59649      unique F087|v59649 -0.924853     FALSE
## F094|v88099        F094      v88099      unique F094|v88099  0.887211      TRUE
## F002|v41283        F002      v41283      unique F002|v41283  0.418866      TRUE
## F028|v17965        F028      v17965      unique F028|v17965  0.562209     FALSE
## F011|v73807        F011      v73807      unique F011|v73807 -0.578273      TRUE
##                   S02       S03       S04       S05       S06       S07
##             <logical> <logical> <logical> <logical> <logical> <logical>
## F048|v14842     FALSE      TRUE     FALSE     FALSE     FALSE     FALSE
## F087|v59649     FALSE     FALSE      TRUE     FALSE     FALSE     FALSE
## F094|v88099     FALSE     FALSE     FALSE     FALSE     FALSE     FALSE
## F002|v41283     FALSE     FALSE     FALSE     FALSE     FALSE     FALSE
## F028|v17965     FALSE      TRUE     FALSE     FALSE     FALSE     FALSE
## F011|v73807     FALSE     FALSE     FALSE     FALSE     FALSE     FALSE
##                   S08       S09       S10 multistateGroup
##             <logical> <logical> <logical>     <character>
## F048|v14842     FALSE     FALSE     FALSE              NA
## F087|v59649     FALSE     FALSE     FALSE              NA
## F094|v88099     FALSE     FALSE     FALSE              NA
## F002|v41283     FALSE     FALSE     FALSE              NA
## F028|v17965     FALSE     FALSE     FALSE              NA
## F011|v73807     FALSE     FALSE     FALSE              NA

Here is a snapshot of the simulation key for QTL simulated as multi-state:

head(rowData(subset(sim, QTL == "multistate")))

## DataFrame with 6 rows and 16 columns
##              feature_id  variant_id         QTL          id mean_beta       S01
##             <character> <character> <character> <character> <numeric> <logical>
## F018|v19728        F018      v19728  multistate F018|v19728 -0.398207      TRUE
## F084|v98536        F084      v98536  multistate F084|v98536  0.524651      TRUE
## F053|v82198        F053      v82198  multistate F053|v82198  0.526131      TRUE
## F060|v42323        F060      v42323  multistate F060|v42323  0.344318     FALSE
## F059|v53070        F059      v53070  multistate F059|v53070 -0.566041      TRUE
## F091|v47866        F091      v47866  multistate F091|v47866 -0.312867     FALSE
##                   S02       S03       S04       S05       S06       S07
##             <logical> <logical> <logical> <logical> <logical> <logical>
## F018|v19728     FALSE      TRUE     FALSE      TRUE      TRUE      TRUE
## F084|v98536     FALSE      TRUE     FALSE      TRUE      TRUE      TRUE
## F053|v82198     FALSE      TRUE     FALSE      TRUE      TRUE      TRUE
## F060|v42323      TRUE     FALSE      TRUE     FALSE     FALSE     FALSE
## F059|v53070     FALSE      TRUE     FALSE      TRUE      TRUE      TRUE
## F091|v47866      TRUE     FALSE      TRUE     FALSE     FALSE     FALSE
##                   S08       S09       S10 multistateGroup
##             <logical> <logical> <logical>     <character>
## F018|v19728     FALSE      TRUE     FALSE          Group1
## F084|v98536     FALSE      TRUE     FALSE          Group1
## F053|v82198     FALSE      TRUE     FALSE          Group1
## F060|v42323      TRUE     FALSE      TRUE          Group2
## F059|v53070     FALSE      TRUE     FALSE          Group1
## F091|v47866      TRUE     FALSE      TRUE          Group2

message("Number of QTL specific to each state-group:")
table(rowData(subset(sim, QTL == "multistate"))$multistateGroup)

## 
## Group1 Group2 
##    218    182

Pairwise sharing

The function plotPairwiseSharing shows the degree of pairwise sharing of significant hits for each combination of two states. Column annotations can be added by specifying a valid column name from the colData of the object.

In the plot below, columns are ordered by the broad cell type (multistateGroup) of the states. We expect states belonging to the same multi-state group to be more related and have a greater degree of sharing of significant eQTLs. Column annotations are used to show the number of significant eQTLs for each state.

sim_sig <- getSignificant(sim)
sim_top <- getTopHits(sim_sig, assay="lfsrs", mode="state")
sim_top <- runPairwiseSharing(sim_top, factor = 0.9)
p1 <- plotPairwiseSharing(sim_top, annotateColsBy=c("nSignificant", "multistateGroup"))
p1

Upset plots

These plots show the set of tests that are significant, but not necessarily shared, by groups of states.

plotUpSet(sim_top, annotateColsBy=c("nSignificant", "multistateGroup"))

Categorizing multi-state QTL tests

Once multi-state test correction is performed, you will want to identify global, multi-state, and unique QTL.

Note that plotCompareStates returns a list with a ggplot2 object as the first element and a table as the second element. These can be accessed using the names “plot” or “counts”.

sim_top <- runTestMetrics(sim_top)

plotCompareStates(sim_top, x="S01", y="S02")

## $plot

## 
## $counts
## 
## both_diverging    both_shared            S01            S02 
##            157            288              5              6

table(rowData(sim_top)$qtl_type)

## 
##     global_diverging        global_shared multistate_diverging 
##                  279                  113                   53 
##    multistate_shared 
##                   11

hist(rowData(sim_top)$nSignificant)

Visualizing multi-state QTL

The function plotQTLClusters can be used to produce a heatmap of the eQTL betas values for each state. Each row is a SNP-gene pair, and columns are states. Row and column annotations can be added by naming column names from the rowData or colData of the input QTLExperiment object.

sim_top_ms <- subset(sim_top, qtl_type_simple == "multistate")

plotQTLClusters(
    sim_top_ms, 
    annotateColsBy=c("multistateGroup"),
    annotateRowsBy=c("qtl_type", "mean_beta", "QTL"))