Using R: plotting the genome on a line
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Imagine you want to make a Manhattan-style plot or anything else where you want a series of intervals laid out on one axis after one another. If it’s actually a Manhattan plot you may have a friendly R package that does it for you, but here is how to cobble the plot together ourselves with ggplot2.
We start by making some fake data. Here, we have three contigs (this could be your chromosomes, your genomic intervals or whatever) divided into one, two and three windows, respectively. Each window has a value that we’ll put on the y-axis.
library(dplyr) library(ggplot2) data <- data_frame(contig = c("a", "a", "a", "b", "b", "c"), start = c(0, 500, 1000, 0, 500, 0), end = c(500, 1000, 1500, 500, 1000, 200), value = c(0.5, 0.2, 0.4, 0.5, 0.3, 0.1))
We will need to know how long each contig is. In this case, if we assume that the windows cover the whole thing, we can get this from the data. If not, say if the windows don’t go up to the end of the chromosome, we will have to get this data from elsewhere (often some genome assembly metadata). This is also where we can decide in what order we want the contigs.
contig_lengths <- summarise(group_by(data, contig), length = max(end))
Now, we need to transform the coordinates on each contig to coordinates on our new axis, where we lay the contings after one another. What we need to do is to add an offset to each point, where the offset is the sum of the lengths of the contigs we’ve layed down before this one. We make a function that takes three arguments: two vectors containing the contig of each point and the position of each point, and also the table of lengths we just made.
flatten_coordinates <- function(contig, coord, contig_lengths) { coord_flat <- coord offset <- 0 for (contig_ix in 1:nrow(contig_lengths)) { on_contig <- contig == contig_lengths$contig[contig_ix] coord_flat[on_contig] <- coord[on_contig] + offset offset <- offset + contig_lengths$length[contig_ix] } coord_flat }
Now, we use this to transform the start and end of each window. We also transform the vector of the length of the contigs, so we can use it to add vertical lines between the contigs.
data$start_flat <- flatten_coordinates(data$contig, data$start, contig_lengths) data$end_flat <- flatten_coordinates(data$contig, data$end, contig_lengths) contig_lengths$length_flat <- flatten_coordinates(contig_lengths$contig, contig_lengths$length, contig_lengths)
It would be nice to label the x-axis with contig names. One way to do this is to take the coordinates we just made for the vertical lines, add a zero, and shift them one position, like so:
axis_coord <- c(0, contig_lengths$length_flat[-nrow(contig_lengths)])
Now it’s time to plot! We add one layer of points for the values on the y-axis, where each point is centered on the middle of the window, followed by a layer of vertical lines at the borders between contigs. Finally, we add our custom x-axis, and also some window dressing.
plot_genome <- ggplot() + geom_point(aes(x = (start_flat + end_flat)/2, y = value), data = data) + geom_vline(aes(xintercept = length_flat), data = contig_lengths) + scale_x_continuous(breaks = axis_coord, labels = contig_lengths$contig, limits = c(0, max(contig_lengths$length_flat))) + xlab("Contig") + ylim(0, 1) + theme_bw()
And this is what we get:
I’m sure your plot will look more impressive, but you get the idea.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.