-
Notifications
You must be signed in to change notification settings - Fork 1
/
genemodel.plot.R
119 lines (98 loc) · 3.87 KB
/
genemodel.plot.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
#'genemodel.plot
#'
#' This function plots a gene model
#'
#' @param model data.frame containing model information. Required columns are "type", "coordinates"
#' @param start start position
#' @param bpstop stop position
#' @param orientation either "foward" or "reverse" indicates the direction of transcription
#' @param xaxis default is TRUE and adds axis above gene model showing position
#' @import graphics
#' @export
#' @examples
#' data(AT5G62640)
#' genemodel.plot(AT5G62640, 25149433, 25152541, "reverse", xaxis=TRUE)
genemodel.plot<-function(model, start, bpstop, orientation, xaxis=TRUE)
{
par(mar=c(1,1,3,1), cex=1)
model<-cbind(model[,1], as.data.frame(stringr::str_split_fixed(model$coordinates, "-", 2)))
colnames(model)<-c("feature", "start", "bpstop")
model$start<-as.numeric(as.character(model$start));model$bpstop<-as.numeric(as.character(model$bpstop))
length<-bpstop-start
if (orientation=="reverse")
{
model$newstart<-bpstop-model$bpstop+1
model$newstop<-bpstop-model$start+1
model<-model[which(model$feature!="exon"),]
model<-model[which(model$feature!="ORF"),]
model<-model[order(model$newstart),]
plot(1, type="l",axes=F,ann=FALSE, xlim=c(start-.03*length, bpstop), ylim=c(-1, .5))
for (i in 2:nrow(model))
{
type<-model$feature[i]
if (type=="coding_region")
{
rect(model$newstart[i], 0, model$newstop[i], .2, col = "steelblue3", border="dodgerblue4", lwd=1)
}
if (type=="intron")
{
mid<-mean(c(model$newstart[i], model$newstop[i]))
segments(x0=model$newstart[i],y0=.1,x1=mid,y1=.2, lwd=1, col="dodgerblue4")
segments(x0=model$newstop[i],y0=.1,x1=mid,y1=.2, lwd=1, col="dodgerblue4")
}
if (type=="3' utr")
{
rect(model$newstart[i], 0, model$newstop[i], .2, col = "lightsteelblue1", border="dodgerblue4", lwd=1)
}
if (type=="5' utr")
{
rect(model$newstart[i], 0, model$newstop[i], .2, col = "lightsteelblue1", border="dodgerblue4", lwd=1)
}
}
x<-c(model$newstop[1], model$newstop[1], model$newstart[1], start-.02*length, model$newstart[1])
y<-c(0,.2,.2,.1,0)
endtype<-model$feature[1]
if (endtype=="coding_region") polygon(x,y, border = "dodgerblue4" , col ="steelblue3" , lwd=1)
else polygon(x,y, border = "dodgerblue4" , col ="lightsteelblue1" , lwd=1)
}
if (orientation=="forward")
{
model$newstart<-start+model$start-1
model$newstop<-start+model$bpstop-1
model<-model[which(model$feature!="exon"),]
model<-model[which(model$feature!="ORF"),]
model<-model[order(model$newstop, decreasing = T),]
plot(1, type="l",yaxt='n',ann=FALSE, xlim=c(start, bpstop+.03*length), ylim=c(-1, .5))
for (i in 2:nrow(model))
{
type<-model$feature[i]
if (type=="coding_region")
{
rect(model$newstart[i], 0, model$newstop[i], .2, col = "steelblue3", border="dodgerblue4", lwd=1)
}
if (type=="intron")
{
mid<-mean(c(model$newstart[i], model$newstop[i]))
segments(x0=model$newstart[i],y0=.1,x1=mid,y1=.2, lwd=1, col="dodgerblue4")
segments(x0=model$newstop[i],y0=.1,x1=mid,y1=.2, lwd=1, col="dodgerblue4")
}
if (type=="3' utr")
{
rect(model$newstart[i], 0, model$newstop[i], .2, col = "lightsteelblue1", border="dodgerblue4", lwd=1)
}
if (type=="5' utr")
{
rect(model$newstart[i], 0, model$newstop[i], .2, col = "lightsteelblue1", border="dodgerblue4", lwd=1)
}
}
x<-c(model$newstart[1], model$newstart[1], model$newstop[1], bpstop+.02*length, model$newstop[1])
y<-c(0,.2,.2,.1,0)
endtype<-model$feature[1]
if (endtype=="coding_region") {polygon(x,y, border = "dodgerblue4" , col ="steelblue3" , lwd=1) }
else {polygon(x,y, border = "dodgerblue4" , col ="lightsteelblue1" , lwd=1)}
}
if (xaxis==T)
{
Axis(side=3, labels=T, cex.axis=0.7)
}
}