#!/usr/bin/Rscript
# File-Name:       prok-genomes-genes-graph.R           
# Date:            2013-01-11
# Author:          James Estevez (User:Estevezj)
# Purpose:         This generates a log-log plot of protein count as a function of genome size.
# Data Used:       ftp://ftp.ncbi.nlm.nih.gov/genomes/GENOME_REPORTS/prokaryotes.txt
# License:         To the extent possible under law, the author(s) have
#                  dedicated all copyright and related and neighboring rights to this software to
#                  the public domain worldwide. This software is distributed without any
#                  warranty. You should have received a copy of the CC0 Public Domain Dedication
#                  along with this software. If not, see
#                  <https://creativecommons.org/publicdomain/zero/1.0/>.
library(grDevices)
library(ggplot2)
library(plyr)
library(taxize)

# Download our tables from NCBI's FTP site. Accessed Fri Jan 11 23:02:49 PST 2013

prok <- read.table("ftp://ftp.ncbi.nlm.nih.gov/genomes/GENOME_REPORTS/prokaryotes.txt", sep="\t", comment.char="!", header=T, stringsAsFactors = F)

prok <- read.table("ncbi-ftp-reports-prokaryotes.txt", sep="\t", comment.char="!", header=T, stringsAsFactors = F)
# Clear missing values ('-')
prok.cut <- prok[(prok$Size..Mb. != '-') & (prok$Proteins != '-'),]

# Set classes
prok.cut$Size..Mb. <- as.numeric(prok.cut$Size..Mb.)
prok.cut$Proteins <- as.numeric(prok.cut$Proteins)
prok.cut$Group <- as.factor(prok.cut$Group)

# From which domain of life does each genome come?
groups <- levels(prok.cut$Group)
get_domain <- function(x){first.hit <- classification(get_uid(x))[[1]] # return the first hit
                          kingdom <- as.character(first.hit[which(first.hit[,"Rank"] == "superkingdom"), 1]) # extract domain
                          return(data.frame(Group = x, Domain = kingdom))
}
domains <- ldply(groups, get_domain)
foo <- prok.cut
prok.cut <- merge(prok.cut, domains, by = "Group")

# Draw our plot
p <- ggplot(prok.cut, aes(Size..Mb., Proteins, color = Domain))

# Save our plot to SVG
svg(filename='ncbi-prok-genomesize.svg', width = 14, height = 8)
p +  geom_point(alpha = 0.5, size = 2) +
  scale_y_log10() +
  scale_x_log10() +
  scale_shape(solid = FALSE) +
  ggtitle("The total genome size and the number of genes in bacteria and archaea.") +
  xlab('Genome size (Megabases)') +
  ylab("Number of protein coding genes") +
  scale_colour_brewer(type="qual", palette=3)
dev.off()