Last updated: 2018-06-06

workflowr checks: (Click a bullet for more information)
  • R Markdown file: up-to-date

    Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

  • Environment: empty

    Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.

  • Seed: set.seed(1)

    The command set.seed(1) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.

  • Session information: recorded

    Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

  • Repository version: 5d5a403

    Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

    Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
    Ignored files:
        Ignored:    .sos/
        Ignored:    data/.sos/
        Ignored:    output/MatrixEQTLSumStats.Portable.Z.coved.K3.P3.lite.single.expanded.V1.loglik.rds
        Ignored:    workflows/.ipynb_checkpoints/
        Ignored:    workflows/.sos/
    Untracked files:
        Untracked:  fastqtl_to_mash_output/
        Untracked:  gtex6_workflow_output/
    Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
Expand here to see past versions:
    File Version Author Date Message
    Rmd b60baa8 Peter Carbonetto 2018-06-06 Misc. revisions to Rmd files.
    Rmd 7c6df92 Peter Carbonetto 2018-06-06 Some misc. updates to the R Markdown files.
    html 542ef1a Peter Carbonetto 2018-06-06 Created new webpage for HeterogeneityTables analysis.
    Rmd 66af923 Peter Carbonetto 2018-06-06 wflow_publish(c(“gtex.Rmd”, “HeterogeneityTables.Rmd”))
    Rmd 9ee901f Peter Carbonetto 2018-06-06 Revised data/results loading steps in HeterogeneityTables.Rmd.
    Rmd 2543833 Peter Carbonetto 2018-06-06 Renamed HeterogeneityTables analysis files.
    html 2543833 Peter Carbonetto 2018-06-06 Renamed HeterogeneityTables analysis files.

Here we summarize overall sharing of effects by sign and by magnitude. Compare the table at the bottom of this page against Table 2 in the manuscript.

Because a major feature of these data is that brain tissues generally show more similar effects than non-brain tissues, we also compute results separately from subsets of brain and non-brain tissues.

Set up environment

First, we load some functions defined for mash analyses.


This is the threshold used to determine which genes have at least one significant effect across tissues.

thresh <- 0.05

Load data and mash results

Load some GTEx summary statistics, as well as some of the results generated from the mash analysis of the GTEx data.

out            <- readRDS("../data/MatrixEQTLSumStats.Portable.Z.rds")
maxbeta        <- out$test.b
maxz           <- out$test.z
standard.error <- out$test.s
out  <- readRDS(paste("../output/MatrixEQTLSumStats.Portable.Z.coved.K3.P3",
                      "lite.single.expanded.V1.posterior.rds",sep = "."))
pm.mash        <- out$posterior.means
pm.mash.beta   <- pm.mash*standard.error
lfsr           <- out$lfsr
lfsr[lfsr < 0] <- 0
tissue.names   <- as.character(read.table("../data/abbreviate.names.txt")[,2])
colnames(lfsr) <- tissue.names

Load the results generated from the mash analysis of the GTEx data after removing the data from the brain tissues.

lfsr.nobrain           <- read.table("../output/nobrainlfsr.txt")[,-1]
colnames(lfsr.nobrain) <- tissue.names[-c(7:16)]
pm.mash.nobrain <-
  as.matrix(read.table("../output/nobrainposterior.means.txt")[,-1]) *

Load the results generated from the mash analysis of the GTEx data for the brain tissues only.

lfsr.brain.only           <- read.table("../output/brainonlylfsr.txt")[,-1]
colnames(lfsr.brain.only) <- tissue.names[c(7:16)]
pm.mash.brain.only <-
  as.matrix(read.table("../output/brainonlyposterior.means.txt")[,-1]) *

Compute overall sharing by sign and magnitude

Compute the amount of eQTL sharing by sign, in all tissues, and separately in brain and non-brain tissues. “Sharing by” sign requires that the effect has the same sign as the strongest effect among tissues.

sigmat  <- (lfsr<=thresh)
nsig    <-  rowSums(sigmat)
signall <- mean(het.norm(pm.mash.beta[nsig>0,])>0)

sigmat          <- (lfsr[,-c(7:16)]<=thresh)
nsig            <-  rowSums(sigmat)
signall.nobrain <- mean(het.norm(pm.mash.beta[nsig,-c(7:16)])>0)

sigmat            <- (lfsr[,c(7:16)]<=thresh)
nsig              <- rowSums(sigmat)
signall.brainonly <- mean(het.norm(pm.mash.beta[nsig>0,c(7:16)])>0)

sigmat      <- (lfsr.nobrain<=thresh)
nsig        <- rowSums(sigmat)
signnobrain <- mean(het.norm(pm.mash.nobrain[nsig>0,])>0)

sigmat        <- (lfsr.brain.only<=thresh)
nsig          <- rowSums(sigmat)
signbrainonly <- mean(het.norm(pm.mash.brain.only[nsig>0,])>0)

Compute the amount of sharing by magnitude, in all tissues, and separately in brain and non-brain tissues. “Sharing by Magnitude” requires that the effect is also within a factor of 2 of the strongest effect.

sigmat <- (lfsr<=thresh)
nsig   <- rowSums(sigmat)
magall <- mean(het.norm(pm.mash.beta[nsig>0,])>0.5)

sigmat <- (lfsr[,-c(7:16)]<=thresh)
nsig   <- rowSums(sigmat)
magall.excludingbrain <- mean(het.norm(pm.mash.beta[nsig>0,-c(7:16)]) > 0.5)

sigmat <- (lfsr[,c(7:16)]<=thresh)
nsig   <- rowSums(sigmat)
magall.brainonly <- mean(het.norm(pm.mash.beta[nsig>0,c(7:16)]) > 0.5)

sigmat     <- (lfsr.nobrain<=thresh)
nsig       <- rowSums(sigmat)
magnobrain <- mean(het.norm(pm.mash.nobrain[nsig>0,]) > 0.5)

sigmat   <- (lfsr.brain.only<=thresh)
nsig     <- rowSums(sigmat)
magbrain <- mean(het.norm(pm.mash.brain.only[nsig>0,]) > 0.5)

Summarize these calculations in a single table. The numbers in parentheses are obtained by the secondary mash analyses on the brain-only and non-brain tissue subsets.

         nrow = 2,ncol = 5,
         dimnames = list(c("shared by sign","shared by magnitude"),
                         c("all tissues","non-brain","(non-brain)",
      digits = 3)
#                     all tissues non-brain (non-brain) brain (brain)
# shared by sign            0.850     0.849       0.882 0.959   0.984
# shared by magnitude       0.359     0.398       0.445 0.764   0.859

The results confirm extensive eQTL sharing among tissues, particularly among the brain tissues; sharing in sign exceeds 85% in all cases, and is as high as 96% among the brain tissues.

Sharing in magnitude is inevitably lower, because sharing in magnitude implies sharing in sign. Overall, on average 36% of tissues show an effect within a factor of 2 of the strongest effect at each top eQTL.

However, within brain tissues this number increases to 76%. That is, not only do eQTLs tend to be shared among the brain tissues, but the effect sizes tend to be quite homogeneous.

Session information

# R version 3.4.3 (2017-11-30)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS High Sierra 10.13.4
# Matrix products: default
# BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
# LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
# locale:
# [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
# attached base packages:
# [1] stats     graphics  grDevices utils     datasets  methods   base     
# loaded via a namespace (and not attached):
#  [1] workflowr_1.0.1.9000 Rcpp_0.12.16         digest_0.6.15       
#  [4] rprojroot_1.3-2      R.methodsS3_1.7.1    backports_1.1.2     
#  [7] git2r_0.21.0         magrittr_1.5         evaluate_0.10.1     
# [10] stringi_1.1.7        whisker_0.3-2        R.oo_1.21.0         
# [13] R.utils_2.6.0        rmarkdown_1.9        tools_3.4.3         
# [16] stringr_1.3.0        yaml_2.1.18          compiler_3.4.3      
# [19] htmltools_0.3.6      knitr_1.20

This reproducible R Markdown analysis was created with workflowr