Comparison of Flaveria trinervia and Flaveria robusta

setwd('~/hydrogen/flaveria/Fr/Fr_2')
fr <- read.csv('Fr_express.tpm', sep="\t", head=T)
key <- read.csv('key_agi.txt',sep='\t', head=F)
names(key) <- c("annotation", "description")

setwd('~/hydrogen/flaveria/Ft/Ft_2')
ft <- read.csv('Ft_express.tpm', sep="\t", head=T)

setwd('~/hydrogen/flaveria/round_robin')
#anno <- read.csv('robin_11.txt', sep="\t", head=F)
anno <- read.csv('flaveria_round_robin_annotation.txt', sep="\t", head=F)
names(anno) <- c("species", "contig", "reference", "annotation", "bitscore", "cascades")
anno$species <- gsub(anno$species, pattern="_cd_hit", replacement="")

library(plyr)
library(reshape2)
library(ggplot2)
fr_anno <- anno[which(anno$species=="Fr"),]
ft_anno <- anno[which(anno$species=="Ft"),]

fr_2 <- merge(fr, fr_anno, by.x="contig", by.y="contig")       # merge tpm data and annotation from rrobin
fr_2 <- fr_2[,-c(2, 21,22)]
key_fr <- fr_2[which(fr_2$annotation %in% key$annotation),]    # select only data that is key
fr_tpm <- ddply(key_fr, .(annotation), numcolwise(sum))        # sum columns, grouping by annotation

ft_2 <- merge(ft, ft_anno, by.x="contig", by.y="contig")
ft_2 <- ft_2[,-c(2, 21,22)]
key_ft <- ft_2[which(ft_2$annotation %in% key$annotation),]
ft_tpm <- ddply(key_ft, .(annotation), numcolwise(sum))        # sum columns, grouping by annotation

fr_tpm <- fr_tpm[,-c(20,21)] # remove unneeded columns from annotation (cascade, bitscore)
ft_tpm <- ft_tpm[,-c(20,21)]

fr_melt <- melt(fr_tpm)

## Using annotation as id variables

ft_melt <- melt(ft_tpm)

## Using annotation as id variables

fr_melt$variable <- gsub(fr_melt$variable, pattern="F[a-z]_[0-9].", replacement="")
ft_melt$variable <- gsub(ft_melt$variable, pattern="F[a-z]_[0-9].", replacement="")

fr_sum <- ddply(fr_melt, .(annotation,variable), numcolwise(mean)) # sum 3 replicate tpms
ft_sum <- ddply(ft_melt, .(annotation,variable), numcolwise(mean))

both <- merge(fr_sum, ft_sum, by=c("annotation","variable"))
names(both) <- c("annotation","section","Fr","Ft")

setwd("~/hydrogen/flaveria")
agis <- key[which(key$annotation %in% both$annotation),]

for (i in 1:nrow(agis)) { 
  tmp <- both[which(both$annotation==paste(agis[i,"annotation"])),]
  c<-ggplot(tmp, aes(x = Fr, y = Ft, group=annotation)) +
     geom_point(size = 2) +
     geom_path() +
     geom_text(aes(label=section)) +
     coord_fixed(ratio = 1, xlim=c(0,max(tmp$Fr,tmp$Ft)), ylim=c(0,max(tmp$Fr,tmp$Ft))) +
     ggtitle(paste("Trinervia Robusta comparison of ", agis[i,"annotation"], " ", agis[i,"description"]))
   print(c)
}