|
0
|
1 #! /usr/bin/env Rscript
|
|
|
2
|
|
|
3 ## FUNCTIONS:
|
|
|
4 letterA <- function(x.pos,y.pos,ht,wt,id=NULL){
|
|
|
5
|
|
|
6 x <- c(0,4,6,10,8,6.8,3.2,2,0,3.6,5,6.4,3.6)
|
|
|
7 y <- c(0,10,10,0,0,3,3,0,0,4,7.5,4,4)
|
|
|
8 x <- 0.1*x
|
|
|
9 y <- 0.1*y
|
|
|
10
|
|
|
11 x <- x.pos + wt*x
|
|
|
12 y <- y.pos + ht*y
|
|
|
13
|
|
|
14 if (is.null(id)){
|
|
|
15 id <- c(rep(1,9),rep(2,4))
|
|
|
16 }else{
|
|
|
17 id <- c(rep(id,9),rep(id+1,4))
|
|
|
18 }
|
|
|
19
|
|
|
20 fill <- c("green","white")
|
|
|
21
|
|
|
22 list(x=x,y=y,id=id,fill=fill)
|
|
|
23 }
|
|
|
24
|
|
|
25 ## T
|
|
|
26 letterT <- function(x.pos,y.pos,ht,wt,id=NULL){
|
|
|
27
|
|
|
28 x <- c(0,10,10,6,6,4,4,0)
|
|
|
29 y <- c(10,10,9,9,0,0,9,9)
|
|
|
30 x <- 0.1*x
|
|
|
31 y <- 0.1*y
|
|
|
32
|
|
|
33 x <- x.pos + wt*x
|
|
|
34 y <- y.pos + ht*y
|
|
|
35
|
|
|
36 if (is.null(id)){
|
|
|
37 id <- rep(1,8)
|
|
|
38 }else{
|
|
|
39 id <- rep(id,8)
|
|
|
40 }
|
|
|
41
|
|
|
42 fill <- "red"
|
|
|
43
|
|
|
44 list(x=x,y=y,id=id,fill=fill)
|
|
|
45 }
|
|
|
46
|
|
|
47 ## C
|
|
|
48 letterC <- function(x.pos,y.pos,ht,wt,id=NULL){
|
|
|
49 angle1 <- seq(0.3+pi/2,pi,length=100)
|
|
|
50 angle2 <- seq(pi,1.5*pi,length=100)
|
|
|
51 x.l1 <- 0.5 + 0.5*sin(angle1)
|
|
|
52 y.l1 <- 0.5 + 0.5*cos(angle1)
|
|
|
53 x.l2 <- 0.5 + 0.5*sin(angle2)
|
|
|
54 y.l2 <- 0.5 + 0.5*cos(angle2)
|
|
|
55
|
|
|
56 x.l <- c(x.l1,x.l2)
|
|
|
57 y.l <- c(y.l1,y.l2)
|
|
|
58
|
|
|
59 x <- c(x.l,rev(x.l))
|
|
|
60 y <- c(y.l,1-rev(y.l))
|
|
|
61
|
|
|
62 x.i1 <- 0.5 +0.35*sin(angle1)
|
|
|
63 y.i1 <- 0.5 +0.35*cos(angle1)
|
|
|
64 x.i1 <- x.i1[y.i1<=max(y.l1)]
|
|
|
65 y.i1 <- y.i1[y.i1<=max(y.l1)]
|
|
|
66 y.i1[1] <- max(y.l1)
|
|
|
67
|
|
|
68 x.i2 <- 0.5 +0.35*sin(angle2)
|
|
|
69 y.i2 <- 0.5 +0.35*cos(angle2)
|
|
|
70
|
|
|
71 x.i <- c(x.i1,x.i2)
|
|
|
72 y.i <- c(y.i1,y.i2)
|
|
|
73
|
|
|
74 x1 <- c(x.i,rev(x.i))
|
|
|
75 y1 <- c(y.i,1-rev(y.i))
|
|
|
76
|
|
|
77 x <- c(x,rev(x1))
|
|
|
78 y <- c(y,rev(y1))
|
|
|
79
|
|
|
80 x <- x.pos + wt*x
|
|
|
81 y <- y.pos + ht*y
|
|
|
82
|
|
|
83 if (is.null(id)){
|
|
|
84 id <- rep(1,length(x))
|
|
|
85 }else{
|
|
|
86 id <- rep(id,length(x))
|
|
|
87 }
|
|
|
88
|
|
|
89 fill <- "blue"
|
|
|
90
|
|
|
91 list(x=x,y=y,id=id,fill=fill)
|
|
|
92 }
|
|
|
93
|
|
|
94
|
|
|
95 ## G
|
|
|
96 letterG <- function(x.pos,y.pos,ht,wt,id=NULL){
|
|
|
97 angle1 <- seq(0.3+pi/2,pi,length=100)
|
|
|
98 angle2 <- seq(pi,1.5*pi,length=100)
|
|
|
99 x.l1 <- 0.5 + 0.5*sin(angle1)
|
|
|
100 y.l1 <- 0.5 + 0.5*cos(angle1)
|
|
|
101 x.l2 <- 0.5 + 0.5*sin(angle2)
|
|
|
102 y.l2 <- 0.5 + 0.5*cos(angle2)
|
|
|
103
|
|
|
104 x.l <- c(x.l1,x.l2)
|
|
|
105 y.l <- c(y.l1,y.l2)
|
|
|
106
|
|
|
107 x <- c(x.l,rev(x.l))
|
|
|
108 y <- c(y.l,1-rev(y.l))
|
|
|
109
|
|
|
110 x.i1 <- 0.5 +0.35*sin(angle1)
|
|
|
111 y.i1 <- 0.5 +0.35*cos(angle1)
|
|
|
112 x.i1 <- x.i1[y.i1<=max(y.l1)]
|
|
|
113 y.i1 <- y.i1[y.i1<=max(y.l1)]
|
|
|
114 y.i1[1] <- max(y.l1)
|
|
|
115
|
|
|
116 x.i2 <- 0.5 +0.35*sin(angle2)
|
|
|
117 y.i2 <- 0.5 +0.35*cos(angle2)
|
|
|
118
|
|
|
119 x.i <- c(x.i1,x.i2)
|
|
|
120 y.i <- c(y.i1,y.i2)
|
|
|
121
|
|
|
122 x1 <- c(x.i,rev(x.i))
|
|
|
123 y1 <- c(y.i,1-rev(y.i))
|
|
|
124
|
|
|
125 x <- c(x,rev(x1))
|
|
|
126 y <- c(y,rev(y1))
|
|
|
127
|
|
|
128 h1 <- max(y.l1)
|
|
|
129 r1 <- max(x.l1)
|
|
|
130
|
|
|
131 h1 <- 0.4
|
|
|
132 x.add <- c(r1,0.5,0.5,r1-0.2,r1-0.2,r1,r1)
|
|
|
133 y.add <- c(h1,h1,h1-0.1,h1-0.1,0,0,h1)
|
|
|
134
|
|
|
135
|
|
|
136
|
|
|
137 if (is.null(id)){
|
|
|
138 id <- c(rep(1,length(x)),rep(2,length(x.add)))
|
|
|
139 }else{
|
|
|
140 id <- c(rep(id,length(x)),rep(id+1,length(x.add)))
|
|
|
141 }
|
|
|
142
|
|
|
143 x <- c(rev(x),x.add)
|
|
|
144 y <- c(rev(y),y.add)
|
|
|
145
|
|
|
146 x <- x.pos + wt*x
|
|
|
147 y <- y.pos + ht*y
|
|
|
148
|
|
|
149
|
|
|
150 fill <- c("orange","orange")
|
|
|
151
|
|
|
152 list(x=x,y=y,id=id,fill=fill)
|
|
|
153
|
|
|
154 }
|
|
|
155
|
|
|
156 Letter <- function(which,x.pos,y.pos,ht,wt){
|
|
|
157
|
|
|
158 if (which == "A"){
|
|
|
159 letter <- letterA(x.pos,y.pos,ht,wt)
|
|
|
160 }else if (which == "C"){
|
|
|
161 letter <- letterC(x.pos,y.pos,ht,wt)
|
|
|
162 }else if (which == "G"){
|
|
|
163 letter <- letterG(x.pos,y.pos,ht,wt)
|
|
|
164 }else if (which == "T"){
|
|
|
165 letter <- letterT(x.pos,y.pos,ht,wt)
|
|
|
166 }else{
|
|
|
167 stop("which must be one of A,C,G,T")
|
|
|
168 }
|
|
|
169
|
|
|
170 letter
|
|
|
171 }
|
|
|
172
|
|
|
173
|
|
|
174
|
|
|
175
|
|
|
176 plot_multiline_logo = function(cons.logo,read=NULL, W=50, setpar=TRUE, gaps = NULL){
|
|
|
177 ## logo - base order - A C G T
|
|
|
178 if (ncol(cons.logo)==5){
|
|
|
179 gaps_prob = cons.logo[,5]
|
|
|
180 }else{
|
|
|
181 gaps_prob = NULL
|
|
|
182 }
|
|
|
183 ps=10 # Point_Size
|
|
|
184 tm=4
|
|
|
185 pwm=as.matrix(cons.logo[,1:4])
|
|
|
186 N=nrow(pwm)
|
|
|
187 Nori=N
|
|
|
188 if (N<W){
|
|
|
189 W=N
|
|
|
190 }
|
|
|
191 s1=seq(1,N,by=W)
|
|
|
192 s2=seq(W,N,by=W)
|
|
|
193 if (length(s2)<length(s1)){
|
|
|
194 pwm=rbind(pwm,matrix(0,nrow=W*length(s1)-N,ncol=4,dimnames=list(NULL,c('A','C','G','T'))))
|
|
|
195 if (!is.null(read)){
|
|
|
196 pwm_read = rbind(read,matrix(0,nrow=W*length(s1)-N,ncol=4,dimnames=list(NULL,c('A','C','G','T'))))
|
|
|
197 }
|
|
|
198 N=nrow(pwm)
|
|
|
199 s2=seq(W,N,by=W)
|
|
|
200 }
|
|
|
201 if (setpar){
|
|
|
202 par(mfrow = c(ceiling(N/W),1), mar=c(1,4,1,0))
|
|
|
203 }
|
|
|
204 for (i in seq_along(s1)){
|
|
|
205 if (!is.null(read)){
|
|
|
206 plot.logo(pwm_read[s1[i]:s2[i],],maxh=2)
|
|
|
207 }
|
|
|
208 plot.logo(pwm[s1[i]:s2[i],],maxh=max(rowSums(cons.logo)))
|
|
|
209 if(!is.null(gaps)){
|
|
|
210 ## transparent rectangles
|
|
|
211 rect((gaps[ ,'start']-s1[i]+1),0, (gaps[,'end']-s1[i]+2), max(pwm), col="#00000005")
|
|
|
212
|
|
|
213 }
|
|
|
214 if(!is.null(gaps_prob)){
|
|
|
215 rect(seq_along(s1[i]:s2[i]),
|
|
|
216 max(rowSums(cons.logo)),
|
|
|
217 seq_along(s1[i]:s2[i])+1,
|
|
|
218 max(rowSums(cons.logo)) - gaps_prob[s1[i]:s2[i]],
|
|
|
219 col="#00000030")
|
|
|
220
|
|
|
221
|
|
|
222 }
|
|
|
223 ticks=intersect(intersect(pretty(pretty(s1[i]:s2[i])+1),s1[i]:s2[i]),1:Nori)
|
|
|
224 axis(1,at=ticks+1.5-s1[i],label=ticks,tick=FALSE)
|
|
|
225 y=pretty(c(0,max(pwm)),n=tm)
|
|
|
226 axis(2,at=y,label=y,las=2,cex.axis=.7)
|
|
|
227 }
|
|
|
228 }
|
|
|
229
|
|
|
230 plot.logo=function(pwm,maxh=NULL){
|
|
|
231 acgt=c("A","C","G","T")
|
|
|
232 pwm = pwm[,acgt]
|
|
|
233 nbp=dim(pwm)[1]
|
|
|
234 if (is.null(maxh)) {maxh=max(rowSums(pwm))}
|
|
|
235
|
|
|
236 plot(0,0,xlim=c(0,nbp),ylim=c(0,maxh),type="n",axes=F,xlab="",ylab="")
|
|
|
237 for ( i in 1:nbp){
|
|
|
238 S=order(pwm[i,])
|
|
|
239 hgts=pwm[i,S]
|
|
|
240 nts=acgt[S]
|
|
|
241 ypos=c(0,cumsum(hgts)[1:3])
|
|
|
242 for (j in 1:4){
|
|
|
243 if (hgts[j]==0) next
|
|
|
244 L=Letter(which=nts[j],x.pos=i,y.pos=ypos[j],ht=hgts[j],wt=1)
|
|
|
245 Id=L$id==1
|
|
|
246 polygon(L$x[Id],L$y[Id],lty=0,col=L$fill[1])
|
|
|
247 if (sum(L$id==2)>0) {
|
|
|
248 polygon(L$x[!Id],L$y[!Id],lty=0,col=L$fill[2])
|
|
|
249 }
|
|
|
250 }
|
|
|
251 }
|
|
|
252 }
|
|
|
253
|
|
|
254
|
|
|
255
|
