-
Notifications
You must be signed in to change notification settings - Fork 11
/
Tspecific.Rmd
94 lines (75 loc) · 2.99 KB
/
Tspecific.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
---
title: "Number of tissue-specific eQTLs in each tissue"
author: "Sarah Urbut, Gao Wang, Peter Carbonetto and Matthew Stephens"
output: workflowr::wflow_html
---
Despite high average levels of sharing of eQTLs among tissues, mash
also identifies eQTLs that are relatively "tissue-specific". Here we
count the number of "tissue-specific" eQTLs in each tissue.
Compare the bar chart at the bottom of this page against Supplementary
Figure 5 in the manuscript.
```{r knitr, message = FALSE, echo = FALSE}
knitr::opts_chunk$set(collapse = TRUE, fig.width = 8,
fig.height = 6,fig.align = "center",
comment = "#")
```
## Set up environment
First, we load some functions defined for mash analyses.
```{r load-functions}
source("../code/normfuncs.R")
```
This is the threshold used to determine which genes have at least one
significant effect across tissues.
```{r 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")
maxb <- 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
lfsr.mash <- out$lfsr
pm.mash.beta <- pm.mash * standard.error
```
For the bar chart below, we use the colours that are
conventionally used to represent the GTEx tissues in plots.
```{r load-tissue-colors}
missing.tissues <- c(7,8,19,20,24,25,31,34,37)
color.gtex <- read.table("../data/GTExColors.txt",sep = '\t',
comment.char = '')[-missing.tissues,]
col = as.character(color.gtex[,2])
```
## Compute number of "tissue-specific" effects for each tissue.
We define "tissue-specific" to mean that the effect is
at least 2-fold larger in one tissue than in any other tissue.
```{r count-tissue-specific-effects}
nsig <- rowSums(lfsr.mash < thresh)
pm.mash.beta.norm <- het.norm(effectsize = pm.mash.beta)
pm.mash.beta.norm <- pm.mash.beta.norm[nsig > 0,]
lfsr.mash <- as.matrix(lfsr.mash[nsig > 0,])
colnames(lfsr.mash) <- colnames(maxz)
a <- which(rowSums(pm.mash.beta.norm > 0.5) == 1)
lfsr.fold <- as.matrix(lfsr.mash[a,])
pm <- as.matrix(pm.mash.beta.norm[a,])
tspec <- NULL
for(i in 1:ncol(pm))
tspec[i] <- sum(pm[,i] > 0.5)
tspec <- as.matrix(tspec)
rownames(tspec) <- colnames(maxz)
```
## Plot number of "tissue-specific" effects for each tissue
```{r plot-counts}
par(mfrow = c(2,1))
barplot(as.numeric(t(tspec)),las = 2,cex.names = 0.75,col = col,
names = colnames(lfsr.fold))
```
Testis stands out as the tissue with the most tissue-specific
effects. Other tissues showing stronger-than-average tissue
specificity include skeletal muscle, thyroid and transformed cell
lines (fibroblasts and LCLs).