-
Notifications
You must be signed in to change notification settings - Fork 11
/
SharingMag.Rmd
78 lines (64 loc) · 2.56 KB
/
SharingMag.Rmd
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
---
title: "Pairwise sharing by magnitude of eQTLs among tissues"
author: "Sarah Urbut, Gao Wang, Peter Carbonetto and Matthew Stephens"
output: workflowr::wflow_html
---
The plot generated here summarizes eQTL sharing by magnitude between
all pairs of tissues. Compare against Figure 6 of the paper.
```{r knitr, message = FALSE, echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, fig.width = 8,
fig.height = 4,fig.align = "center",
comment = "#")
```
## Set up environment
First, we load the lattice package used for generating the plot below.
```{r load-pkgs, message = FALSE}
library(lattice)
```
## Load data and mash results
In the next code chunk, we load some GTEx summary statistics, as well
as some of the results generated from the mash analysis of the GTEx
data.
```{r load-results}
out <- readRDS("../data/MatrixEQTLSumStats.Portable.Z.rds")
maxb <- out$test.b
maxz <- out$test.z
out <-readRDS(paste("../output/MatrixEQTLSumStats.Portable.Z.coved.K3.P3",
"lite.single.expanded.V1.posterior.rds",sep = "."))
pm.mash <- out$posterior.means
lfsr.all <- out$lfsr
standard.error <- maxb/maxz
pm.mash.beta <- pm.mash*standard.error
```
## Compute sharing-by-magnitude statistics
For every pair of tissues, we count the proportion of effects
significant in either tissue that are within 2-fold magnitude of one
another.
```{r generate-matrix}
thresh <- 0.05
pm.mash.beta <- pm.mash.beta[rowSums(lfsr.all<0.05)>0,]
lfsr.mash <- lfsr.all[rowSums(lfsr.all<0.05)>0,]
shared.fold.size <- matrix(NA,nrow = ncol(lfsr.mash),ncol=ncol(lfsr.mash))
colnames(shared.fold.size) <- rownames(shared.fold.size) <- colnames(maxz)
for (i in 1:ncol(lfsr.mash))
for (j in 1:ncol(lfsr.mash)) {
sig.row=which(lfsr.mash[,i]<thresh)
sig.col=which(lfsr.mash[,j]<thresh)
a=(union(sig.row,sig.col))
quotient=(pm.mash.beta[a,i]/pm.mash.beta[a,j])
shared.fold.size[i,j] = mean(quotient > 0.5 & quotient < 2)
}
```
## Plot heatmap of sharing by magnitude
Generate the heatmap using the "levelplot" function from the lattice
package.
```{r heatmap-sharing-magnitude, fig.height=10, fig.width=10}
all.tissue.order <- read.table("../data/alltissueorder.txt")[,1]
clrs <- colorRampPalette(rev(c("#D73027","#FC8D59","#FEE090","#FFFFBF",
"#E0F3F8","#91BFDB","#4575B4")))(64)
lat <- shared.fold.size[rev(all.tissue.order),rev(all.tissue.order)]
lat[lower.tri(lat)] <- NA
n <- nrow(lat)
print(levelplot(lat[n:1,],col.regions = clrs,xlab = "",ylab = "",
colorkey = TRUE))
```