-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy patheda_chrom-states.Rmd
78 lines (52 loc) · 1.53 KB
/
eda_chrom-states.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
---
title: "eda_chrom-states.Rmd"
author: "Jason Torres"
date: "April 25, 2017"
output: html_document
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
# Setup
```{r}
# Approach originally proposed by Parker et al. 2013. PNAS
"%&%" <- function(a,b) paste0(a,b)
library("dplyr")
library("data.table")
pre <- "/Users/jtorres/FUSE/"
in.dir <- pre %&% "reference/islet/"
cur.dir <- in.dir %&% "stretch_enhancers/"
cs.df <- fread(in.dir %&% "Pancreat_islet_15_dense.reformatted_colours.bed",sep="\t")
cs.df$len <- cs.df$V3 - cs.df$V4
```
Write subsetted file with only weak,strong enhancer states (6,7,8,9)
```{r}
sub.df <- filter(cs.df,V4==6 | V4==7 | V4==8 | V4==9) %>% select(one_of("V1","V2","V3","V4"))
write.table(sub.df,file=cur.dir%&%"Pancreat_islet_15_dense.en-only.bed",col.names=F,row.names=F,
sep="\t",quote=F)
```
Create merged enhancer file
```{bash}
sh merge.sh
```
Get distribution
```{r}
en.df <- fread(cur.dir %&% "islet_merged_enhancers.bed")
en.df$len <- en.df$V3 - en.df$V2
hist(en.df$len,breaks=100)
summary(en.df$len/1000)
qua <- quantile(en.df$len/1000,probs=seq(0,1,0.05))
q90 <- qua[grepl(pattern="90%",names(qua))] # 4.2 kb
st.df <- filter(en.df,len>=q90*1000)
write.table(st.df,file=cur.dir%&%"Pancreat_islet_15_dense.en-only.perc90.bed",col.names=F,row.names=F,
sep="\t",quote=F)
```
Make stretch enhancer bed
```{bash}
sh merge2.sh
```
Evalute stretch enhancer bed
```{r}
s.df <- fread(cur.dir %&% "Pancreat_stretch-enhancers.bed")
s.df$len <- s.df$V3 - s.df$V2
```