# plot the results of splice site prediction
#define some colours
rgb <- c("#009E73", "#D55E00", "#0072B2")
# make two graphs on top of each other
par(mfrow = c(2, 1))
# First consider the 5' splice site prediction and read the
# output from
# the Perl code
data <- read.table("score5.txt", sep = "\t", header = TRUE)
# for the plot, we need to know about the sequence length
seqlen <- max(data$pos)
# for the plot we need to know about the maximum score
max_score <- max(data$score)
# make a plot for the 5' splice site data
plot(0, type = "n", lwd = 2, xlim = c(0, seqlen),
ylim = c(0, max_score * 1.1), main = "Splice site scoring",
xlab = "Position", ylab = "Score")
#print a legend
legend(seqlen * 0.7, max_score, "5prime", col = rgb[2],
lwd = 1)
# plot the splice site scores
for (i in (1:seqlen)) {
lines(c(data$pos[i], data$pos[i]), c(0, data$score[i]), col = rgb[2],
lw = 2)
}
# plot the location of exons (we do not know these from the
# prediction)
lines(c(268, 331), c(max_score/2, max_score/2), col = rgb[1],
lw = 4)
lines(c(447, 1054), c(max_score/2, max_score/2), col = rgb[1],
lw = 4)
# Now consider the 3' splice site and read the output from
# the Perl code
data <- read.table("score3.txt", sep = "\t", header = TRUE)
max_score <- max(data$score)
plot(0, type = "n", lwd = 2, xlim = c(0, seqlen),
ylim = c(0, max_score * 1.1), xlab = "Position", ylab = "Score")
legend(seqlen * 0.7, max_score, "3prime", col = rgb[3],
lwd = 1)
for (i in (1:seqlen)) {
lines(c(data$pos[i], data$pos[i]), c(0, data$score[i]), col = rgb[3],
lw = 2)
}
# plot the location of exons as in the previous graph
lines(c(268, 331), c(max_score/2, max_score/2), col = rgb[1],
lw = 4)
lines(c(447, 1054), c(max_score/2, max_score/2), col = rgb[1],
lw = 4)