## Viricel A, Pante E, Dabin W, Simon-Bouhet B. Applicability of RAD-tag genotyping for inter-familial comparisons: empirical data from two cetaceans ## Note that this script was made for a Unix system ## This script was used to make histograms of the Blast2GO annotation results. ############################################################ ## get data setwd("~/Dropbox/macmito/temp_amelia/blast2go_with_xml/") read.csv("All_annotation_results_08-13-2013.csv", header=T) -> ann head(ann) ############################################################ ## pivot tables keep = 10 droplevels(subset(ann, ann$GO.Group == "P")) -> subann # Biological Process table(subann$Term, subann$Set) -> term.per.pop.P droplevels(subset(ann, ann$GO.Group == "F")) -> subann # Molecular Function table(subann$Term, subann$Set) -> term.per.pop.F droplevels(subset(ann, ann$GO.Group == "C")) -> subann # Cellular Component table(subann$Term, subann$Set) -> term.per.pop.C head(term.per.pop.C, keep) ; head(term.per.pop.F, keep) ; head(term.per.pop.P, keep) ############################################################ ## Venn Diagrams: library(limma) vennCounts(term.per.pop.P) -> venn.P vennCounts(term.per.pop.F) -> venn.F vennCounts(term.per.pop.C) -> venn.C par(mfrow=c(1,3)) ; cx=1.2 vennDiagram(venn.P, cex=cx, main="Biological Process") vennDiagram(venn.F, cex=cx, main="Molecular Function") vennDiagram(venn.C, cex=cx, main="Cellular Component") ############################################################ ## Bar Graphs go = 1000# number of GO categories to display term.per.pop.P[order(term.per.pop.P[,1], decreasing=T),1][1:go] -> PcCons # count, Biological Process, conserved as.vector(term.per.pop.P[order(term.per.pop.P[,2], decreasing=T),2][1:go]) -> PcVar term.per.pop.F[order(term.per.pop.F[,1], decreasing=T),1][1:go] -> FcCons # count, Molecular Function, conserved as.vector(term.per.pop.F[order(term.per.pop.F[,2], decreasing=T),2][1:go]) -> FcVar term.per.pop.C[order(term.per.pop.C[,1], decreasing=T),1][1:go] -> CcCons # count, Cellular Component, conserved as.vector(term.per.pop.C[order(term.per.pop.C[,2], decreasing=T),2][1:go]) -> CcVar txt=0.7 par(las=1, mai=c(0.35,3,0.2,0.1), mfrow=c(3,2)) barplot(PcCons, cex.names=txt, main="Biological Process, Conserved tags", horiz=TRUE) barplot(PcVar, cex.names =txt, main="Biological Process, Variable tags", horiz=TRUE) barplot(FcCons, cex.names =txt, main="Molecular Function, Conserved tags", horiz=TRUE) barplot(FcVar, cex.names =txt, main="Molecular Function, Variable tags", horiz=TRUE) barplot(CcCons, cex.names =txt, main="Cellular Component, Conserved tags", horiz=TRUE) barplot(CcVar, cex.names =txt, main="Cellular Component, Variable tags", horiz=TRUE)