Comparison of genes with tissue-specific eQTLs against other genes

Load data and mash results

In the next code chunk, we load some GTEx summary statistics (average gene expression values and z-scores), as well as some of the results generated from the mash analysis of the GTEx data.

Expression is here defined as median across individuals of the log Reads per Kilobase Mapped (RPKM).

data <- read.csv("../data/AvgExpr.csv.gz",header = TRUE)
out  <- readRDS("../data/MatrixEQTLSumStats.Portable.Z.rds")
maxz <- out$strong.z
qtl.names <-
  sapply(1:length(rownames(maxz)),
         function(x) unlist(strsplit(rownames(maxz)[x], "[_]"))[[1]])
rownames(data) <- data[,1]
expr.data      <- data[,-1]
out    <- readRDS(paste("../output/MatrixEQTLSumStats.Portable.Z.coved.K3.P3",
                        "lite.single.expanded.V1.posterior.rds",sep = "."))
pmash <- out$posterior.means
lfsr  <- out$lfsr
colnames(lfsr) <- colnames(pmash) <- colnames(maxz)
rownames(lfsr) <- rownames(pmash) <- rownames(maxz)
expr.sort <- expr.data[rownames(expr.data)%in%qtl.names,]
a         <- match(qtl.names,rownames(expr.sort))
expr.sort <- expr.sort[a,]
missing.tissues    <- c(7,8,19,20,24,25,31,34,37)
exp.sort           <- expr.sort[,-missing.tissues]
colnames(exp.sort) <- colnames(maxz)

Define plotting functions

Here we define a couple functions used to compare the densities and CDFs of gene expression levels.

plot_tissuespecifictwo = function(tissuename,lfsr,ymax,curvedata,title,
                                  thresh=0.05,subset=1:44,exclude=0.01) {
  index_tissue=which(colnames(lfsr) %in% tissuename);
  ybar=title
  
  # Create a matrix showing whether or not lfsr satisfies threshold.
  sigmat = lfsr <= thresh;
  sigs=which(rowSums(sigmat[,index_tissue,drop=FALSE])==length(tissuename) &
             rowSums(sigmat[,-index_tissue,drop=FALSE])==0)
  norm.spec=density(curvedata[sigs,index_tissue])
  norm.other=density(curvedata[-sigs,index_tissue])
  xmin=min(norm.spec$x,norm.other$x)-1
  ymin=min(norm.spec$y,norm.other$y)
  xmax=max(norm.spec$x,norm.other$x)+1

  plot(norm.other,xlim = c(xmin,xmax),ylim=c(0,ymax),
       col="black",main=tissuename,xlab="log(RPKM)")
  
  lines(norm.spec,col="red")
  legend("right",legend = c("All Genes","Tissue Specific"),
         col=c("black","red"),pch=c(1,1))
  }

plot_tissuespecificthree = function(tissuename,lfsr,ymax,curvedata,title,
                                    thresh=0.05,subset=1:44,exclude=0.01) {
  index_tissue=which(colnames(lfsr) %in% tissuename);
  ybar=title
  
  # Create a matrix showing whether or not lfsr satisfies threshold.
  sigmat = lfsr <= thresh
  sigs=which(rowSums(sigmat[,index_tissue,drop=FALSE])==length(tissuename) &
             rowSums(sigmat[,-index_tissue,drop=FALSE])==0)
  norm.spec=ecdf(curvedata[sigs,index_tissue])
  norm.other=ecdf(curvedata[-sigs,index_tissue])
  plot(norm.other,ylim=c(0,ymax),
       col="black",main=tissuename,xlab="log(RPKM)")
  
    lines(norm.spec,col="red",cex=0.1)
    legend("right",legend = c("All Genes","Tissue Specific"),
           col=c("black","red"),pch=c(1,1))
  } 

Distribution of expression levels in testis

The two plots below compare the densities and cumulative distribution functions of the expression levels for all genes (black), and for genes identified as having a “tissue-specific” eQTL (red) in testis.

plot_tissuespecifictwo(tissuename = "Testis",lfsr = lfsr,
                       curvedata = log(exp.sort),title = "Test",
                       thresh = 0.05 ,ymax=0.5)

plot_tissuespecificthree(tissuename = "Testis",lfsr = lfsr,
                         curvedata = log(exp.sort),title = "Test",
                         thresh = 0.05 ,ymax=1)

The distribution functions are reasonably similar.

Distribution of expression levels in thyroid

Next we show the same two plots for thyroid.

plot_tissuespecifictwo(tissuename = "Thyroid",lfsr = lfsr,
                       curvedata = log(exp.sort),title = "Test",
                       thresh = 0.05,ymax = 0.5)

plot_tissuespecificthree(tissuename = "Thyroid",lfsr = lfsr,
                         curvedata = log(exp.sort),title = "Test",
                         thresh = 0.05,ymax = 1)

Distribution of expression levels in whole blood

Next, we look at the distributions in whole blood cells.

plot_tissuespecifictwo(tissuename = "Whole_Blood",lfsr = lfsr,
                       curvedata = log(exp.sort),title = "Test",
                       thresh = 0.05,ymax=0.5)

plot_tissuespecificthree(tissuename = "Whole_Blood",lfsr = lfsr,
                         curvedata = log(exp.sort),title = "Test",
                         thresh = 0.05 ,ymax=1)

Distribution of expression levels in fibroblasts

Finally, we examine the gene expression distributions in fibroblasts.

plot_tissuespecifictwo(tissuename = "Cells_Transformed_fibroblasts",
                       lfsr = lfsr,curvedata = log(exp.sort),title = "Test",
                       thresh = 0.05,ymax=0.5)

plot_tissuespecificthree(tissuename = "Cells_Transformed_fibroblasts",
                        lfsr = lfsr,curvedata = log(exp.sort),title = "Test",
                        thresh = 0.05 ,ymax=1)

In each case, the distribution functions are similar. This tells us that tissue-specific eQTLs are not simply reflecting tissue-specific expression.

Session information

sessionInfo()
# R version 3.4.3 (2017-11-30)
# Platform: x86_64-apple-darwin15.6.0 (64-bit)
# Running under: macOS High Sierra 10.13.5
# 
# 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.17         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