Comparison of Flaveria bidentis 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/Fb/Fb_2')
fb <- read.csv('Fb_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"),]
fb_anno <- anno[which(anno$species=="Fb"),]

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

fb_2 <- merge(fb, fb_anno, by.x="contig", by.y="contig")
fb_2 <- fb_2[,-c(2, 21,22)]
key_fb <- fb_2[which(fb_2$annotation %in% key$annotation),]
fb_tpm <- ddply(key_fb, .(annotation), numcolwise(sum))        # sum columns, grouping by annotation

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

fr_melt <- melt(fr_tpm)

## Using annotation as id variables

fb_melt <- melt(fb_tpm)

## Using annotation as id variables

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

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

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

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 = Fb, group=annotation)) +
     geom_point(size = 2) +
     geom_path() +
     geom_text(aes(label=section)) +
     coord_fixed(ratio = 1, xlim=c(0,max(tmp$Fr,tmp$Fb)), ylim=c(0,max(tmp$Fr,tmp$Fb))) +
     ggtitle(paste("Bidentis Robusta comparison of ", agis[i,"annotation"], " ", agis[i,"description"]))
   print(c)
}