iCLIP/RBP_occupancy_norm.R at master · xxy103/iCLIP · GitHub

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
library(RColorBrewer)
library(scales)

col_set <- brewer.pal(n = 8, 'Set1')

smooth_occu <- function(x, step=3) {
  y <- x
  x <- c(x[1:step], x)
  for (i in 1:length(y)) {
    y[i] <- mean(x[i:(i+step)])
  }
  return(y)
}

ko_si_5ss <- read.table('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_silenced_downstream_5SS_KO_dist2occu')
wt_si_5ss <- read.table('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_silenced_downstream_5SS_WT_dist2occu')
ko_en_5ss <- read.table('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_enhanced_downstream_5SS_KO_dist2occu')
wt_en_5ss <- read.table('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_enhanced_downstream_5SS_WT_dist2occu')
ko_si_3ss <- read.table('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_silenced_upstream_3SS_KO_dist2occu')
wt_si_3ss <- read.table('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_silenced_upstream_3SS_WT_dist2occu')
ko_en_3ss <- read.table('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_enhanced_upstream_3SS_KO_dist2occu')
wt_en_3ss <- read.table('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_enhanced_upstream_3SS_WT_dist2occu')


#pdf('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_v2.pdf')
pdf('occupancy/MeCP2_KO_RBFOX2_KI_SE_RPKM_2_v6.pdf', hei=5)
par(mfrow=c(2, 2))
#par(mar=c(5,4,4,1))
par(mar=c(2,4,1,1))
# enhanced SE events
plot(wt_en_3ss$V1, smooth_occu(wt_en_3ss$V2, step=10), ylim=c(0,1), type='l', bty='l', las=2, lwd=1.5, col=col_set[1], ylab='Normalized crosslink events', xlab='', xaxt='n')
lines(ko_en_3ss$V1, smooth_occu(ko_en_3ss$V2, step=10), col=col_set[2], lwd=1.5)
axis(1, at=c(0,50,100,150), lab=c(-100, -50, '3SS', 50))
legend('topleft', c('WT', 'KO'), col = col_set[1:2], lwd = 3, bty='n')
abline(h=0)
#rect(100, -0.02, 200, 0.02, col='gray', border = NA)
rect(100, -0.03, 200, 0.03, col='black', border = NA)
x <- smooth_occu(wt_en_3ss$V2, step=10)
y <- smooth_occu(ko_en_3ss$V2, step=10)
x <- c(x, x[-(1:120)])
y <- c(y, y[-(1:120)])
p <- c()
fc <- c()
for (i in 1:150) {
  p[i] <- t.test(x[i:(i+30)], y[i:(i+30)])$p.value
  fc[i] <- log2(mean(x[i:(i+30)])/mean(y[i:(i+30)]))
}
start <- which(p.adjust(p) < 0.01 & fc > log2(1.5))
end <- which(p.adjust(p) < 0.01 & fc > log2(1.5))+30
rect(min(start), 0, max(end), 10, col = alpha(col=col_set[1], 0.1), border = NA)

#par(mar=c(5,1,4,4))
par(mar=c(2,1,1,4))
plot(wt_en_5ss$V1, smooth_occu(wt_en_5ss$V2, step=10), ylim=c(0,1), type='l', bty='n', yaxt='n', las=2, lwd=1.5, col=col_set[1], ylab='', xlab='', main='', xaxt='n')
lines(ko_en_5ss$V1, smooth_occu(ko_en_5ss$V2, step=10), col=col_set[2], lwd=1.5)
axis(1, at=c(-100,0,50,100,150,200), lab=c('', -50, '5SS', 50, 100, ''))
abline(h=0)
#rect(-50, -0.02, 50, 0.02, col='gray', border = NA)
rect(-50, -0.03, 50, 0.03, col='black', border = NA)
x <- smooth_occu(wt_en_5ss$V2, step=10)
y <- smooth_occu(ko_en_5ss$V2, step=10)
x <- c(x, x[-(1:120)])
y <- c(y, y[-(1:120)])
p <- c()
fc <- c()
for (i in 1:150) {
  p[i] <- t.test(x[i:(i+30)], y[i:(i+30)])$p.value
  fc[i] <- log2(mean(x[i:(i+30)])/mean(y[i:(i+30)]))
}
start <- which(p.adjust(p) < 0.01 & fc > log2(1.5))
end <- which(p.adjust(p) < 0.01 & fc > log2(1.5))+30
rect(min(start), 0, max(end), 10, col = alpha(col=col_set[1], 0.1), border = NA)

#par(mar=c(5,4,4,1))
par(mar=c(2,4,1,1))
# silenced SE events
plot(wt_si_3ss$V1, smooth_occu(wt_si_3ss$V2, step=10), ylim=c(0,1), type='l', bty='l', las=2, lwd=1.5, col=col_set[3], ylab='Normalized crosslink events', xlab='', xaxt='n')
lines(ko_si_3ss$V1, smooth_occu(ko_si_3ss$V2, step=10), col=col_set[4], lwd=1.5)
axis(1, at=c(0,50,100,150), lab=c(-100, -50, '3SS', 50))
legend('topleft', c('WT', 'KO'), col = col_set[3:4], lwd = 3, bty='n')
abline(h=0)
#rect(100, -0.02, 200, 0.02, col='gray', border = NA)
rect(100, -0.03, 200, 0.03, col='black', border = NA)
x <- smooth_occu(wt_si_3ss$V2, step=10)
y <- smooth_occu(ko_si_3ss$V2, step=10)
x <- c(x, x[-(1:120)])
y <- c(y, y[-(1:120)])
p <- c()
fc <- c()
for (i in 1:150) {
  p[i] <- t.test(x[i:(i+30)], y[i:(i+30)])$p.value
  fc[i] <- log2(mean(x[i:(i+30)])/mean(y[i:(i+30)]))
}
start <- which(p.adjust(p) < 0.01 & fc > log2(1.5))
end <- which(p.adjust(p) < 0.01 & fc > log2(1.5))+30
rect(min(start), 0, max(end), 10, col = alpha(col=col_set[1], 0.1), border = NA)

#par(mar=c(5,1,4,4))
par(mar=c(2,1,1,4))
plot(wt_si_5ss$V1, smooth_occu(wt_si_5ss$V2, step=10), ylim=c(0,1), type='l', bty='n', yaxt='n', las=2, lwd=1.5, col=col_set[3], ylab='', xlab='', main='', xaxt='n')
lines(ko_si_5ss$V1, smooth_occu(ko_si_5ss$V2, step=10), col=col_set[4], lwd=1.5)
axis(1, at=c(-100,0,50,100,150,200), lab=c('', -50, '5SS', 50, 100, ''))
abline(h=0)
#rect(-50, -0.02, 50, 0.02, col='gray', border = NA)
rect(-50, -0.03, 50, 0.03, col='black', border = NA)
x <- smooth_occu(wt_si_5ss$V2, step=10)
y <- smooth_occu(ko_si_5ss$V2, step=10)
x <- c(x, x[-(1:120)])
y <- c(y, y[-(1:120)])
p <- c()
fc <- c()
for (i in 1:150) {
  p[i] <- t.test(x[i:(i+30)], y[i:(i+30)])$p.value
  fc[i] <- log2(mean(x[i:(i+30)])/mean(y[i:(i+30)]))
}
start <- which(p.adjust(p) < 0.01 & fc > log2(1.5))
end <- which(p.adjust(p) < 0.01 & fc > log2(1.5))+30
rect(min(start), 0, max(end), 10, col = alpha(col=col_set[1], 0.1), border = NA)

dev.off()