-
Notifications
You must be signed in to change notification settings - Fork 11
/
SharingHist.Rmd
112 lines (92 loc) · 3.32 KB
/
SharingHist.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
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
---
title: "Add appropriate title here."
author: "Sarah Urbut"
output: workflowr::wflow_html
---
*Add text here.*
```{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 some functions defined for the mash analyses.
```{r load-functions}
source("../code/normfuncs.R")
```
*Add text here.*
```{r define-thresh}
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.
```{r load-results}
out <- readRDS("../data/MatrixEQTLSumStats.Portable.Z.rds")
maxbeta <- out$test.b
maxz <- out$test.z
out <- readRDS(paste("../output/MatrixEQTLSumStats.Portable.Z.coved.K3.P3",
"lite.single.expanded.V1.posterior.rds",sep = "."))
lfsr <- out$lfsr
pm.mash <- out$posterior.means
standard.error <- maxbeta/maxz
pm.mash.beta <- pm.mash * standard.error
```
# Add section title here
Let's plot heterogneity by magnitude from the global analysis.
```{r hindex-plot-magnitude}
sigmat <- (lfsr<=thresh)
nsig <- rowSums(sigmat)
par(mar = c(4,4,2,1))
par(oma = c(8,4,0,0) + 0.1)
par(mfrow = c(1,3))
hist((het.func(het.norm(effectsize=pm.mash.beta[nsig>0,]),threshold=0.5)),
main="",xlab="",breaks=0.5:44.5,col="grey",freq=FALSE,ylim=c(0,0.9),
xaxt="n")
axis(1,at = seq(1, 44, by=1),labels = c(1:44))
mtext("All Tissues")
sigmat <- (lfsr[,-c(7:16)] <= thresh)
nsig <- rowSums(sigmat)
hist((het.func(het.norm(effectsize=pm.mash.beta[nsig>0,-c(7:16)]),
threshold=0.5)),main="",xlab="",breaks=0.5:34.5,col="grey",
freq=FALSE,ylim=c(0,0.9),xaxt="n")
axis(1,at = seq(1, 34, by=1),labels = c(1:34))
mtext("Non-Brain Tissues")
sigmat <- (lfsr[,c(7:16)]<=thresh)
nsig <- rowSums(sigmat)
brain.norm <- het.norm(effectsize=pm.mash.beta[nsig>0,c(7:16)])
hist(het.func(brain.norm,threshold=0.5),main="",xlab="",breaks=0.5:10.5,
col="grey",freq=FALSE,xaxt="n",ylim=c(0,0.9))
axis(1, at=seq(1, 10, by=1), labels=c(1:10))
mtext("Brain Tissues")
````
## Add section title here
Now, let's make the same plot with all tissue effects measuring the
number of tissues which have a sign equivalent to max effect:
````{r hindex-plot-sign}
sign.func <- function (normeffectsize)
apply(normeffectsize,1,function(x)(sum(x>0)))
sigmat <- (lfsr<=thresh)
nsig <- rowSums(sigmat)
par(mar = c(4,4.5,2,1))
par(oma = c(8,4,0,0) + 0.1)
par(mfrow = c(1,3))
hist(sign.func(het.norm(effectsize=pm.mash.beta[nsig>0,])),main="",xlab="",
breaks=0.5:44.5,col="grey",freq=FALSE,xaxt="n",ylim=c(0,0.9))
axis(1, at=seq(1, 44, by=1), labels=c(1:44))
mtext("All Tissues")
sigmat <- (lfsr[,-c(7:16)] <= thresh)
nsig <- rowSums(sigmat)
hist(sign.func(het.norm(effectsize = pm.mash.beta[nsig>0,-c(7:16)])),
main="",xlab="",breaks=0.5:34.5,col="grey",freq=FALSE,ylim=c(0,0.9),
xaxt="n")
axis(1, at=seq(1, 34, by=1), labels=c(1:34))
mtext("Non-Brain Tissues")
sigmat <- (lfsr[,c(7:16)]<=thresh)
nsig <- rowSums(sigmat)
brain.norm <- het.norm(effectsize=pm.mash.beta[nsig>0,c(7:16)])
hist(sign.func(brain.norm),main="",xlab="",breaks=0.5:10.5,col="grey",
freq=FALSE,xaxt="n",ylim=c(0,0.9))
axis(1, at=seq(1, 10, by=1), labels=c(1:10))
mtext("Brain Tissues")
```