# Download our tables from NCBI's FTP site. Accessed 14:30PST, 18 December 2012
prok <- read.table("ftp://ftp.ncbi.nlm.nih.gov/genomes/GENOME_REPORTS/prokaryotes.txt", sep="\t", comment.char="!", header=T)

# Pull release dates, while dropping rows lacking a release date.
prok  <- as.Date(prok$Release.Date[prok$Release.Date != '-'],format="%Y/%m/%d")

# Bin our dates by month and year, tabulate, and save to a dataframe.
prok.cut <- as.data.frame( 
table(
  as.Date(
    cut(prok, "month")
  )
)
)

# Correct our column titles, calculate a running total, and reconvert from factor to date
colnames(prok.cut) <- c("Date", "Total")
prok.cut$Total <- cumsum(prok.cut$Total)
prok.cut$Date <- as.Date(prok.cut$Date)

# DNA Sequencing Costs from NHGRI: http://www.genome.gov/sequencingcosts/
# Data from http://www.genome.gov/pages/der/sequencing_cost.pptx

# After munging the pptx, download the tables from pastebin. Accessed 12:42PST, 2012-12-20
seq.cost <- read.table("http://pastebin.com/raw.php?i=NA6c4i70", header=TRUE)

# Format the date.
seq.cost$Date  <- as.Date(seq.cost$Date,format="%m-%d-%Y")

# Draw our plots
library("ggplot2")
library("grid")
library("scales")

(p <- ggplot(prok.cut, aes(Date, Total)) +  geom_area() +  ggtitle("Bacterial and archeal genome sequences submitted to Genbank") + xlab('Time') +  ylab("Total number of genomes")
)
(mb <- ggplot(seq.cost, aes(Date, USD.per.Mb)) + geom_point(colour = "blue") +
stat_smooth(color="#984EA3")+  
ggtitle("Cost to sequence one million nucleotides") +
xlab('Time') +
ylab("USD per MB") +
scale_y_continuous(labels = dollar)
)
(genome <- ggplot(seq.cost, aes(Date, USD.per.Genome)) + geom_point(colour = "red") +
stat_smooth(method='lm',color="#FC8D62")+  
ggtitle("Cost to sequence one human genome") +
xlab('Time') +
ylab("USD per genome") +
scale_y_log10(labels = dollar)
)
# This part is based on Hadley's Ggplot2 book (doi:10.1007/978-0-387-98141-3_8)
# Save our plot to SVG
library(grDevices)
svg(filename='ncbi-genomes.svg', width = 15, height = 10)
grid.newpage()
pushViewport(viewport(layout = grid.layout(2, 2)))
vplayout <- function(x, y)
viewport(layout.pos.row = x, layout.pos.col = y)
print(p, vp = vplayout(1, 1:2))
print(mb, vp = vplayout(2, 1))
print(genome, vp = vplayout(2, 2))
dev.off()