Visualize CAMI assemblies on radar chart

additional_analyses/comparison/comparison.R

@@ -302,4 +302,86 @@ dev.off()
 #save.image("/home/shaman/Work/Data/integrated-omics-pipeline/MS_analysis/comparison/data_usage/comparison_revision.Rdat")
 #load("/home/shaman/Work/Data/integrated-omics-pipeline/MS_analysis/comparison/data_usage/comparison_revision.Rdat")
 
+###############################################################################
+## CAMI data analyses
+###############################################################################
+
+indir <- "/home/shaman/Work/Data/integrated-omics-pipeline/MS_analysis/metaquast_output/CAMI"
+samples <- c("CAMI_medium", "CAMI_low")
+
+for(i in seq_along(samples)){
+   dat <- read.delim(paste(indir, "/", samples[i],  "/", "transposed_report.tsv", sep=""), header=T, stringsAsFactors=F) 
+
+   assign(paste(samples[i], "_quast", sep=""), dat )
+   rm(dat)
+}
+
+###############################################################################################################################
+### FUNCTION: 3-axis radar chart
+plot.3axis <- function(dat, cols, dens, font, linetype, 
+			     fsize, linewd, mcex, lwd, plwd, title = "") {
+dat <- CAMI_low_quast
+
+dat <- dat[,c("Assembly", 
+       "X..contigs.....1000.bp.", 
+       "N50", 
+       "X..predicted.genes..unique."
+       )]
+
+colnames(dat)  <- c("Assembly", "contigs \u2265 1kb", 
+		    "N50 length", "no. of genes")
+
+dat$Assembly <- factor(dat$Assembly, levels=c(as.character(dat$Assembly), "max", "min")) 
+dat <- rbind(
+	     c("max", apply(dat[, 2:ncol(dat)], 2, max)), 
+	     c("min", apply(dat[, 2:ncol(dat)], 2, min)), 
+	     dat)
+
+dat[,2:ncol(dat)] <- sapply(dat[,2:ncol(dat)], as.numeric)
+	     
+# Calculate maximum and minimum for positive values
+dat[1, !names(dat) %in% c("Assembly")] <- 
+    dat[1, !names(dat) %in% c("Assembly")] * 1.025
+dat[2, !names(dat) %in% c("Assembly")] <- 
+    dat[2, !names(dat) %in% c("Assembly")] * 0.75
+
+par(mar=c(0,0,8,0)) 
+radarchart(dat[,-1], pcol=cols , pfcol=dens, plty=linetype, plwd=plwd,
+	   cglwd=lwd, palcex=font, calcex=font, vlcex=font, cglcol="darkgray",
+	   paxislabels = c("Assembly", "contigs \u2265 1kb", "N50 length", "no. of genes"))
+	   oma=c(0,0,0,0)#, title = title, cex.main = 6)
+}
+
+cols <- makeTransparent("darkblue", "darkgreen", "gold", alpha=0.75)
+dens <- makeTransparent("blue", "green", "gold",  alpha=0.10)
+font=10
+linetype <- c(5,5,5)
+fsize=10
+linewd=15
+mcex=8
+lwd=16
+plwd=15
+
+
+pdf("/home/shaman/Documents/Publications/IMP-manuscript/figures/third_iteration/CAMI_low_radarChart.pdf", 
+    width=38, height=25)
+plot.3axis(CAMI_low_quast, cols, dens, font, linetype, fsize, linewd, mcex, lwd, plwd) 
+dev.off()
+
+pdf("/home/shaman/Documents/Publications/IMP-manuscript/figures/third_iteration/CAMI_medium_radarChart.pdf", 
+    width=38, height=25)
+plot.3axis(CAMI_medium_quast, cols, dens, font, linetype, fsize, linewd, mcex, lwd, plwd) 
+dev.off()
+
+## Generate legend for figure
+legend_labels <- c("IMP-megahit", "MOCAT","Gold_standard")
+
+pdf("/home/shaman/Documents/Publications/IMP-manuscript/figures/third_iteration/CAMI_legend.pdf", width=45, height=25)
+plot(1, type="n", axes=FALSE, xlab="", ylab="")
+legend("bottom", legend=legend_labels, xpd = TRUE, horiz = TRUE,
+       bty = "o", col = cols, lty=linetype,
+       lwd = c(8,8,8,8), cex = 5, box.lty=1, box.lwd=0,
+       box.col="gray", text.col=cols, text.font=c(2,2,2,2))
+dev.off()
+