-
Notifications
You must be signed in to change notification settings - Fork 1
/
plot_CG.meth.r
185 lines (168 loc) · 4.35 KB
/
plot_CG.meth.r
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
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
#!/usr/bin/env Rscript
options(scipen = 999)
library(ggplot2)
library(tidyverse)
library(dplyr)
library(reshape2)
##### on chr11 #####
files <- list.files(pattern="*.mESC.*chr11.CG.meth.xls")
meth <- data.frame()
for(i in files){
tmp <- read_delim(i, delim="\t")
tmp$smp <- i
meth <- rbind(meth, tmp)
}
meth.w <- spread(meth[,c(1,2,10)],key=smp,value=mods)
p <- ggplot(meth, aes(x = pos, y = mods, color = smp)) +
geom_bar(stat="identity") +
facet_wrap(~smp,ncol=1) +
theme_minimal() +
ylab("meth") +
theme(legend.position="none")
ggsave(
paste0("chr11.meth.pdf"),
p,
width = 15,
height = 8 ,
dpi = 300
)
##### on chr13 #####
files <- list.files(pattern="*.mESC.*chr13.CG.meth.xls")
meth <- data.frame()
for(i in files){
tmp <- read_delim(i, delim="\t")
tmp$smp <- i
meth <- rbind(meth, tmp)
}
meth.w <- spread(meth[,c(1,2,10)],key=smp,value=mods)
p <- ggplot(selmeth, aes(x = pos, y = mods, color = smp)) +
geom_bar(stat="identity") +
facet_wrap(~smp,ncol=1) +
theme_minimal() +
ylab("meth") +
theme(legend.position="none")
ggsave(
paste0("chr13.meth.pdf"),
p,
width = 15,
height = 5 ,
dpi = 300
)
meth.w <- spread(meth[,c(1,2,10)],key=smp,value=mods)
selmeth <- meth.w[,c(1,3,5)] %>% as.data.frame()
p <- ggplot(selmeth, aes(x = selmeth[,2], y = selmeth[,3])) +
geom_point(color="#F0C415") +
theme_minimal() +
ylab("SMRT-TAPS") +
xlab("Nano-TAPS") +
theme(legend.position="none")
ggsave(
paste0("chr13.meth.cor.pdf"),
p,
width = 3,
height = 3 ,
dpi = 300
)
cor(selmeth[,2],selmeth[,3])
##### hbv #####
files <- list.files(pattern="*.hbv.*taps.CG.meth.xls")
meth <- data.frame()
for(i in files){
tmp <- read_delim(i, delim="\t")
tmp$smp <- i
meth <- rbind(meth, tmp)
}
selmeth <- meth[meth$varratio<0.2, ]
meth.w <- spread(selmeth[,c(1,2,10)],key=smp,value=mods)
meth.w$nanopore_huh1.vs.notaps <- meth.w$nanopore_rep3.hbv_huh1.taps.CG.meth.xls - meth.w$nanopore_rep3.hbv_huh1.no_taps.CG.meth.xls
meth.w$pacbio_huh1.vs.notaps <- meth.w$pacbio_rep2.hbv_huh1.taps.CG.meth.xls - meth.w$pacbio_rep2.hbv_huh1.no_taps.CG.meth.xls
selmeth <- meth.w[,c(1,8,9,12,13,16,17)] %>% melt(id.vars=("pos"))
selmeth$value[selmeth$value<0] <- 0
p <- ggplot(selmeth, aes(x = pos, y = value, color = variable)) +
geom_bar(stat="identity") +
facet_wrap(~variable,ncol=1) +
theme_minimal() +
ylab("meth") +
ylim(0,1)+
theme(legend.position="none")
ggsave(
paste0("hbv.meth.huh1_all.pdf"),
p,
width = 15,
height = 12 ,
dpi = 300
)
selmeth <- meth[(meth$fwdmC+1)/(meth$revmC+1)>0.25 & (meth$fwdmC+1)/(meth$revmC+1)<4 & meth$varratio <0.2, ]
p <- ggplot(selmeth, aes(x = pos, y = mods, color = smp)) +
geom_bar(stat="identity") +
facet_wrap(~smp,ncol=1) +
theme_minimal() +
ylab("meth") +
theme(legend.position="none")
ggsave(
paste0("hbv.meth.pdf"),
p,
width = 15,
height = 12 ,
dpi = 300
)
selmeth <- selmeth[grep("hbv_huh1.taps",selmeth$smp),]
p <- ggplot(selmeth, aes(x = pos, y = mods, color = smp)) +
geom_bar(stat="identity") +
facet_wrap(~smp,ncol=1) +
theme_minimal() +
ylab("meth") +
scale_color_manual(values=c("#046699","#EFC319"))+
theme(legend.position="none") +
xlim(0,3248)
ggsave(
paste0("hbv.meth.huh1_taps.pdf"),
p,
width = 15,
height = 12 ,
dpi = 300
)
##### lambda #####
files <- list.files(pattern="*.lambda.*taps.CG.meth.xls")
meth <- data.frame()
for(i in files){
tmp <- read_delim(i, delim="\t")
tmp$smp <- i
meth <- rbind(meth, tmp)
}
meth.w <- spread(meth[,c(1,2,10)],key=smp,value=mods)
p <- ggplot(meth, aes(x = pos, y = mods, color = smp)) +
geom_bar(stat="identity") +
facet_wrap(~smp,ncol=1) +
theme_minimal() +
ylab("meth") +
theme(legend.position="none")
ggsave(
paste0("lambda.meth.pdf"),
p,
width = 15,
height = 5 ,
dpi = 300
)
##### 4kb #####
files <- list.files(pattern="*.4kb.*.CG.meth.xls")
meth <- data.frame()
for(i in files){
tmp <- read_delim(i, delim="\t")
tmp$smp <- i
meth <- rbind(meth, tmp)
}
meth.w <- spread(meth[,c(1,2,10)],key=smp,value=mods)
p <- ggplot(meth, aes(x = pos, y = mods, color = smp)) +
geom_bar(stat="identity") +
facet_wrap(~smp,ncol=1) +
theme_minimal() +
ylab("meth") +
theme(legend.position="none")
ggsave(
paste0("4kb.meth.pdf"),
p,
width = 15,
height = 5 ,
dpi = 300
)