new dginn file format 4 bats + mondrian plots

covid_comp.Rnw

@@ -341,6 +341,21 @@ n12, n23, n13, n123,
 category=c("DGINN-Young-primate", "Young-primate", "Cooper-primate"))
 @
 
+\subsection{Mondrian}
+
+<<mondrianprimates>>=
+
+library(Mondrian)
+
+monddata<-as.data.frame(tmp$Gene.name)
+monddata$primates_dginn_young<-ifelse(tmp$PosSel_PamlM7M8=="Y", 1,0)
+monddata$primates_young<-ifelse(tmp$pVal.M8vsM7<0.05, 1, 0)
+monddata$primates_cooper<-ifelse(tmp$cooper.primates.M7.M8_p_val<0.05, 1, 0)
+
+mondrian(monddata[,2:4])
+@
+
+
 %\subsection{Comparaison des codons?}
 
 %Subtable with lines with both methods showing positive selection.
@@ -405,9 +420,21 @@ listjoined<-mapply(c, listdginn, listjanet, SIMPLIFY=FALSE)
 
 Lecture du tableau.
 
+<<readtab4test, eval=F>>=
+# Tableau tel qu'il est à cet étape du script, pour travailler sur l'ajout du nouveau tableau bats sans recommencer au début à chaque fois.
+tab<-read.delim("COVID_PAMLresults_332hits_plusBatScreens_plusDGINN_20200506.txt", 
+	header=TRUE, sep="\t")
+@
+
 <<>>=
-dginnbats<-read.delim("/home/adminmarie/Documents/CIRI_BIBS_projects/2020_05_Etienne_covid/data/2020_bats_completeResults.csv", 
+#dginnbatsold<-read.delim("/home/adminmarie/Documents/CIRI_BIBS_projects/2020_05_Etienne_covid/data/2020_bats_completeResults.csv", 
+#		  fill=T, h=T)
+
+
+dginnbats<-read.delim("/home/adminmarie/Documents/CIRI_BIBS_projects/2020_05_Etienne_covid/data/DGINN_202005281339summary_cleaned.tab", 
 		  fill=T, h=T)
+
+
 dim(dginnbats)
 names(dginnbats)
 length(unique(dginnbats$Gene))
@@ -425,42 +452,9 @@ Merge tables:
 
 <<bats>>=
 
-add_col<-function(dginnfile=dginnbats, method="PamlM1M2", prefixe="bats_"){
-
-tmp<-dginn[dginnfile$Method==method,
-	   c("Gene", "Omega", "PosSel", "PValue", "NbSites", "PSS")]
-
-names(tmp)<-c("cooper.batsGene", paste0(prefixe, "Omega_", method), 
-	      paste0(prefixe, "PosSel_", method), paste0(prefixe, "PValue_", method), 
-	      paste0(prefixe, "NbSites_", method), paste0(prefixe, "PSS_", method))
-
-tab<-merge(tab, tmp, by="cooper.batsGene", all.x=T)
-
-return(tab)
-}
-
-tab<-add_col(dginnfile=dginnbats, method="PamlM1M2", prefixe="bats_")
+names(dginnbats)<-c("File", "bats_Name", "cooper.batsGene", paste0("bats_", names(dginnbats)[-(1:3)]))
 
-tab<-add_col(dginnfile=dginnbats, method="PamlM7M8", prefixe="bats_")
-
-tab<-add_col(dginnfile=dginnbats, method="BppM1M2", prefixe="bats_")
-
-tab<-add_col(dginnfile=dginnbats, method="BppM7M8", prefixe="bats_")
-
-
-head(na.omit(tab[, c("cooper.batsGene", "cooper.batsAverage_dNdS","cooper.batsM7.M8_p_value", "bats_Omega_PamlM7M8","bats_PValue_PamlM7M8")]))
-
-# Manip pour la colonne BUSTED
-
-tmp<-dginn[dginn$Method=="BUSTED",c("Gene", "Omega", "PosSel", "PValue")]
-names(tmp)<-c("cooper.batsGene", "bats_Omega_BUSTED", "bats_PosSel_BUSTED", "bats_PValue_BUSTED")
-tab<-merge(tab, tmp, all.x=T, by="cooper.batsGene")
-
-tmp<-dginn[dginn$Method=="MEME",c("Gene", "NbSites", "PSS")]
-names(tmp)<-c("cooper.batsGene", "bats_NbSites_MEME", "bats_PSS_MEME")
-tab<-merge(tab, tmp, all.x=T, by="cooper.batsGene")
-
-dim(tab)
+tab<-merge(tab,dginnbats, by="cooper.batsGene", all.x=T)
 @
 
 
@@ -470,22 +464,23 @@ dim(tab)
 
 <<omegaM7M8bats>>=
 
-plot(tab$cooper.batsAverage_dNdS, tab$bats_Omega_PamlM7M8, 
+plot(tab$cooper.batsAverage_dNdS, tab$bats_omegaM0codeml, 
      xlab="Omega Cooper-bats", ylab="Omega DGINN-bats")
 @
 
 \subsubsection{pvalues pour M7M8}
 
 <<pvalM7M8bats>>=
+tab$bats_codemlM7M8.p.value<-as.numeric(as.character(tab$bats_codemlM7M8.p.value))
 
-plot(tab$cooper.batsM7.M8_p_value, tab$bats_PValue_PamlM7M8, pch=20, 
+plot(tab$cooper.batsM7.M8_p_value, tab$bats_codemlM7M8.p.value, pch=20, 
      xlab="p-value Cooper-bats", ylab="p-value DGINN-bats", main="M7vM8 Paml")
 
-points(tab$cooper.batsM7.M8_p_value[tab$cooper.batsM7.M8_p_value>0.05 & tab$bats_PValue_PamlM7M8<0.05], 
-       tab$bats_PValue_PamlM7M8[tab$cooper.batsM7.M8_p_value>0.05 & tab$bats_PValue_PamlM7M8<0.05], 
+points(tab$cooper.batsM7.M8_p_value[tab$cooper.batsM7.M8_p_value>0.05 & tab$bats_codemlM7M8.p.value<0.05], 
+       tab$bats_codemlM7M8.p.value[tab$cooper.batsM7.M8_p_value>0.05 & tab$bats_codemlM7M8.p.value<0.05], 
        col="red", pch=20)
-points(tab$cooper.batsM7.M8_p_value[tab$cooper.batsM7.M8_p_value<0.05 & tab$bats_PValue_PamlM7M8>0.05], 
-       tab$bats_PValue_PamlM7M8[tab$cooper.batsM7.M8_p_value<0.05 & tab$bats_PValue_PamlM7M8>0.05], 
+points(tab$cooper.batsM7.M8_p_value[tab$cooper.batsM7.M8_p_value<0.05 & tab$bats_codemlM7M8.p.value>0.05], 
+       tab$bats_codemlM7M8.p.value[tab$cooper.batsM7.M8_p_value<0.05 & tab$bats_codemlM7M8.p.value>0.05], 
        col="green", pch=20)      
 
 
@@ -531,17 +526,17 @@ legend("topleft", c("<0.05 in Hawkins but >0.05 in Cooper",
 
 <<pvalM7M8compautre>>=
 
-plot(tab$hawkins_Positive.Selection..M8vM8a.p.value, tab$bats_PValue_PamlM7M8, 
+plot(tab$hawkins_Positive.Selection..M8vM8a.p.value, tab$bats_codemlM7M8.p.value, 
 	pch=20, xlab="p-value hawkins-bats", ylab="p-value DGINN-bats", main="M7vM8 Paml")
 
 points(tab$hawkins_Positive.Selection..M8vM8a.p.value[tab$hawkins_Positive.Selection..M8vM8a.p.value>0.05 &
-	tab$bats_PValue_PamlM7M8<0.05], 
-    tab$bats_PValue_PamlM7M8[tab$hawkins_Positive.Selection..M8vM8a.p.value>0.05 &
-    tab$bats_PValue_PamlM7M8<0.05], col="red", pch=20)
+	tab$bats_codemlM7M8.p.value<0.05], 
+    tab$bats_codemlM7M8.p.value[tab$hawkins_Positive.Selection..M8vM8a.p.value>0.05 &
+    tab$bats_codemlM7M8.p.value<0.05], col="red", pch=20)
 points(tab$hawkins_Positive.Selection..M8vM8a.p.value[tab$hawkins_Positive.Selection..M8vM8a.p.value<0.05 & 
-	tab$bats_PValue_PamlM7M8>0.05], 
-    tab$bats_PValue_PamlM7M8[tab$hawkins_Positive.Selection..M8vM8a.p.value<0.05 & 
-    tab$bats_PValue_PamlM7M8>0.05], col="green", pch=20)      
+	tab$bats_codemlM7M8.p.value>0.05], 
+    tab$bats_codemlM7M8.p.value[tab$hawkins_Positive.Selection..M8vM8a.p.value<0.05 & 
+    tab$bats_codemlM7M8.p.value>0.05], col="green", pch=20)      
 
 
 legend("topleft", c("<0.05 in DGINN-bats but >0.05 in Hawkins",
@@ -561,7 +556,7 @@ I will draw a venn diagramm for the positive genes in the 3 analyses.
 \subsubsection{subtab}
 
 <<subbats>>=
-tmp<-na.omit(tab[,c("Gene.name", "bats_PValue_PamlM7M8", "hawkins_Positive.Selection..M8vM8a.p.value", "cooper.batsM7.M8_p_value")])
+tmp<-na.omit(tab[,c("Gene.name", "bats_codemlM7M8.p.value", "hawkins_Positive.Selection..M8vM8a.p.value", "cooper.batsM7.M8_p_value")])
 dim(tmp)
 @
 
@@ -569,24 +564,43 @@ dim(tmp)
 
 \subsubsection{figure}                                               
 <<vennbats>>=
-area1dginn<-sum(tmp$bats_PValue_PamlM7M8<0.05, na.rm=T)
+area1dginn<-sum(tmp$bats_codemlM7M8.p.value<0.05, na.rm=T)
 area2hawk<-sum(tmp$hawkins_Positive.Selection..M8vM8a.p.value<0.05, na.rm=T)
 area3coop<-sum(tmp$cooper.batsM7.M8_p_value<0.05, na.rm=T)
 
 
-n12<-sum(tmp$bats_PValue_PamlM7M8<0.05 & tmp$hawkins_Positive.Selection..M8vM8a.p.value<0.05, na.rm=T)
+n12<-sum(tmp$bats_codemlM7M8.p.value<0.05 & tmp$hawkins_Positive.Selection..M8vM8a.p.value<0.05, na.rm=T)
 
 n23<-sum(tmp$hawkins_Positive.Selection..M8vM8a.p.value<0.05 & tmp$cooper.batsM7.M8_p_value<0.05, na.rm=T)
 
-n13<-sum(tmp$bats_PValue_PamlM7M8<0.05 & tmp$cooper.batsM7.M8_p_value<0.05, na.rm=T)
+n13<-sum(tmp$bats_codemlM7M8.p.value<0.05 & tmp$cooper.batsM7.M8_p_value<0.05, na.rm=T)
+
+
+n123<-sum(tmp$bats_codemlM7M8.p.value<0.05 & tmp$hawkins_Positive.Selection..M8vM8a.p.value<0.05 & tmp$cooper.batsM7.M8_p_value<0.05, na.rm=T)
+
+draw.triple.venn(area1dginn, area2hawk, area3coop, 
+n12, n23, n13, n123, 
+category=c("DGINN-Young-bats", "Hawkins-bats", "Cooper-bats"))
+
 
+@
 
-n123<-sum(tmp$bats_PValue_PamlM7M8<0.05 & tmp$hawkins_Positive.Selection..M8vM8a.p.value<0.05 & tmp$cooper.batsM7.M8_p_value<0.05, na.rm=T)
+\subsection{Mondrian}
 
+<<mondrianbats>>=
+library(Mondrian)
 
-draw.triple.venn(area1dginn, area2hawk, area3coop, n12, n23, n13, n123, category=c("dginn", "hawkins", "cooper"))
+monddata<-as.data.frame(tmp$Gene.name)
+monddata$bats_dginn<-ifelse(tmp$bats_codemlM7M8.p.value<0.05, 1,0)
+monddata$bats_hawkins<-ifelse(tmp$hawkins_Positive.Selection..M8vM8a.p.value<0.05, 1, 0)
+monddata$bats_cooper<-ifelse(tmp$cooper.batsM7.M8_p_value<0.05, 1, 0)
+
+mondrian(monddata[,2:4])
 @
 
+\section{To do}
+
+Comparaison G4 pas G4
 
 \end{document}