How can we plot alignment distribution ?

Question

[MAPK]

I am looking at figure 8G in this article http://journals.plos.org/plospathogens/article?id=10.1371/journal.ppat.1005640 .I was wondering if anyone has idea on how to get alignment distribution plotted (perhaps using R)? Thanks for your help.

Answer 1

[Kevin Blighe]

They don't give much away in the methods, just stating: "A density map of reads along the sense and antisense strands of the virus was created in R." - they do not mention any specific R package. The figure legend itself gives more information:

Analysis of small RNA read counts along the viral genome, as a function of RNA-seq levels. Small RNA RPM counts per nucleotide were determined by strand then divided by the average RPM counts for the RNA-seq for that genotype and strand. Graphical representation of the virus is in the middle, with the two ORFs indicated, and genome coordinates are along the top-most edge. Sense strand reads are plotted in blue above the genome figure and antisense reads are plotted in orange below. The scale for sense strand and antisense strand values is shown as density heatmaps above the plot.

Looking at it, I would be able to generate that using base R functions. Each rectangle would be just a different elongated plot without axis labels or margins and plotted side by side using the mfrow parameter supplied to par (). The coloured lines would be vertical ablines whose colour is shaded by intensity of reads/expression.

Here is a fairly comprehensive starting point fr you. I believe that someone in the group of Sean Davis also developed a R package that does something similar, and there are undoubtedly many other ways too. To adapt it, you'll have to scale the positioning of the vertical ablines to reflect the viral genome, and you'll have to use colorRampPalette to create a colour gradient for each abline based on reads / expression.

pdf("test.pdf", width=9, height=6)
  par(mar=c(0,4,0,0), mfrow=c(7,1))

  require(RColorBrewer)
  require(plotrix)
  plot.new()
  color.legend(0.1, 0.45, 0.4, 0.65, cex=1.2, legend=c("-1","0","1"),   rect.col=colorRampPalette(c("darkblue","blue2","blue","white","red","red2","darkred"))(1000))
  color.legend(0.6, 0.45, 0.9, 0.65, cex=1.2, legend=c("-1","0","1"),  rect.col=colorRampPalette(c("darkblue","blue2","blue","white","red","red2","darkred"))(1000))

  par(xpd=TRUE)

  plot.new() ; abline(v=seq(0,1, 0.2), lwd=4, col="red") ; box(lwd=2, col="grey75") ; mtext(side=2, cex=1.0, "Sample\nX", font=2)
  plot.new() ; abline(v=seq(0,1, 0.25), lwd=4, col="firebrick1") ; box(lwd=2, col="grey75") ; mtext(side=2, cex=1.0, "Sample\nY", font=2)
  plot.new() ; abline(v=seq(0,1, 0.3), lwd=4, col="skyblue") ; box(lwd=2, col="grey75") ; mtext(side=2, cex=1.0, "Sample\nZ", font=2)
  plot.new() ; abline(v=seq(0,1, 0.1), lwd=4, col="gold") ; box(lwd=2, col="grey75") ; mtext(side=2, cex=1.0, "Sample\nZ", font=2)
  plot.new() ; abline(v=seq(0,1, 0.4), lwd=4, col="limegreen") ; box(lwd=2, col="grey75") ; mtext(side=2, cex=1.0, "Sample\nControl 1", font=2)
  plot.new() ; abline(v=seq(0,1, 0.05), lwd=4, col="royalblue") ; box(lwd=2, col="grey75") ; mtext(side=2, cex=1.0, "Sample\nControl 2", font=2)
dev.off()