-
Notifications
You must be signed in to change notification settings - Fork 3
/
5abcd.Rmd
91 lines (74 loc) · 3.67 KB
/
5abcd.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
---
title: "Figure 5"
output:
html_document:
df_print: paged
---
```{r eval=TRUE, message=FALSE, warning=FALSE, include=TRUE}
suppressPackageStartupMessages(source("scripts/themes.R"))
```
```{r eval=FALSE, message=FALSE, warning=FALSE, include=TRUE}
tm.facs <- loadObject(filename = url("https://scfind.cog.sanger.ac.uk/indexes/tm_facs.rds"))
atac <- loadObject(filename = url("https://scfind.cog.sanger.ac.uk/indexes/atacseq.rds"))
```
## Figure 5a - This will generate .bed format file for the UCSC Genome Browser
```{r eval=FALSE, message=FALSE, warning=FALSE, include=TRUE}
source("scripts/6ac.R")
for( tissue in c("Heart", "Liver", "Spleen"))
{
generateUCSCBedFileForEnhancers(atac, query.gene = "Nr4a1",
tissue = tissue, window = c(1e5, 0),
target.window = c("chr15", 101097277, 101105225 ),
annotation = loadAnnotation("data/fig6ac.tsv"))
}
```
## Figure 5b - This will generate .bed format file for the UCSC Genome Browser
```{r eval=FALSE, message=FALSE, warning=FALSE, include=TRUE}
for( tissue in c("Heart", "Liver", "Spleen"))
{
generateUCSCBedFileForEnhancers(atac, query.gene = "Cd36",
tissue = tissue, window = c(1e5, 0),
target.window = c("chr5", 17188577, 17388756 ),
annotation = loadAnnotation("data/fig6ac.tsv"))
}
```
## Figure 5c
```{r, warning=FALSE}
df.6d <- read.table("data/fig6d.csv", sep = ",", header = T, stringsAsFactors = F)
p.6d <- ggplot(df.6d, aes(x=celltype, y=motif, fill=density)) +
theme_minimal() + theme(axis.title.y = element_blank(), axis.title.x = element_blank(), axis.text.x=element_text(hjust = 0, angle = -30), text = element_text(size = 14), strip.text.x = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
labs(fill = "Number of\nmotif/peak") +
scale_fill_gradient2(high = colorRampPalette(rev(brewer.pal(n = 7, name ="RdPu")))(1)) +
facet_wrap(~ tissue, scales = 'free_x') + geom_tile() +
geom_text(aes(label = pval))
p.6d_ <- ggarrange(
p.6d + theme(legend.position = "none"),
ggarrange(get_legend(p.6d),
bottom = text_grob("p-value\n*\t<= 0.05\n**\t<= 0.01\n***\t<= 0.001"),
ncol = 1, nrow = 3),
ncol = 2, widths = c(3,1)
)
```
## Figure 5d
```{r, warning=FALSE}
df <- readRDS("data/fig5d.rds")
p.5d <- ggplot(df, aes( x = dist_tss, y = fraction, colour = spec_celltype)) + geom_point(size = 1) + facet_wrap(~ spec_celltype, ncol = 1) +
theme_minimal() + scale_x_log10(breaks = trans_breaks("log10", function(x) 10^x), labels = trans_format("log10", math_format(10^.x))) + labs(x = "Distance from TSS (bp)", y = "Fraction of cells") +
theme( panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
panel.background = element_blank(),
legend.position="none",
axis.text = element_text(size = 18),
axis.text.x = element_text( size = 18 ),
axis.title = element_text( size = 18) ,
strip.text.x = element_text(size = 18)
) + scale_colour_manual(values = c(cbPalette, tritPalette, RColorBrewer::brewer.pal(8,colorcode))) +
geom_text_repel(data=subset(df,dist_tss >= 5e5 | fraction >= .15), aes(dist_tss, fraction, label = key , colour = spec_celltype),color = "black", size = 5, box.padding = unit(2, "lines")
)
```
```{r echo=FALSE, fig.height=7, fig.retina=10, fig.width=7, message=FALSE, warning=FALSE}
p.6d_
```
```{r echo=FALSE, fig.height=8, fig.retina=10, fig.width=7, message=FALSE, warning=FALSE}
p.5d
```