diff --git a/R/assessCallingPerformance.R b/R/assessCallingPerformance.R
index a2b49604f..bddfcb8ab 100644
--- a/R/assessCallingPerformance.R
+++ b/R/assessCallingPerformance.R
@@ -1,4 +1,16 @@
-#d <- read.table("sim_calls.table", header=T)
+require("plotrix")
+args = commandArgs(TRUE);
+
+onCMDLine = ! is.na(args[1])
+
+file = "sim_calls.table"
+info = "interactive R"
+if ( onCMDLine ) {
+    file = args[1]
+    d <- read.table(file, header=T)
+    pdf(args[2])
+    info = args[3]
+}
 
 d$sim.VAR <- d$sim.AC > 0
 d$called.VAR <- d$called.AC > 0
@@ -8,13 +20,23 @@ MODES = unique(d$sim.MODE)
 NS = unique(d$called.AN / 2)
 DEPTHS = unique(d$sim.DP)
 
+addSection <- function(name) {
+    par("mar", c(5, 4, 4, 2))
+    frame()
+    title(name, cex=2)
+}   
+
+addSection(paste("Calling performance report: nSamples = ", NS, "\n  info:", info))
+
 results <- expand.grid(Q = QS, mode = MODES, nSamples = NS, depth = DEPTHS)
 results$sensitivity = 0
 results$specificity = 0
 
 determineRates <- function(raw, Q, mode, depth) {
     sub <- subset(raw, sim.Q == Q & sim.MODE == mode & sim.DP == depth)
-    ct <- table(sub$called.VAR, sub$sim.VAR, dnn = c("called.VAR", "sim.VAR"))
+    print(c(Q,mode,depth, dim(sub)))
+    ct <- table(sub$called.VAR, sub$sim.VAR, dnn = c("called.VAR", "sim.VAR"), useNA = "always")
+    print(ct)
     sensitivity = ct[2,2] / sum(ct[,2]) 
     specificity = ct[1,1] / sum(ct[,1])
     list(sensitivity = sensitivity, specificity = specificity, ct = ct)
@@ -28,7 +50,7 @@ for ( i in 1:(dim(results)[1]) ) {
 }
 
 for ( depth in DEPTHS ) 
-    boxplot(called.AC ~ sim.AC, data = subset(d, called.DP == depth * NS), main = paste("Depth of coverage ", depth), xlab = "Simulation AC", ylab = "Called AC")
+    boxplot(called.AC ~ sim.AC, data = subset(d, called.DP == depth * NS), main = paste("Depth of coverage ", depth), xlab = "Simulation AC", ylab = "Called AC", outwex=0.5, col = "cornflowerblue")
 
 print(results)
 
@@ -36,6 +58,14 @@ par(mfcol=c(2,1))
 for ( Qt in QS ) {
     x <- subset(results, Q == Qt)
     print(x)
-    plot(x$depth, x$sensitivity, type="b")
-    plot(x$depth, x$specificity, type="b")
+    plot(x$depth, x$sensitivity, type="b", main = paste("Q score", Qt), xlab = "Depth", ylab="Sensitivity")
+    plot(x$depth, x$specificity, type="b", xlab = "Depth", ylab="Specificity")
 }
+
+par(mfcol=c(1,1))
+plot(0,0, type="n", frame.plot=F, ann=F, axes=F)
+addtable2plot(-1, -1, data.frame(Q=results$Q, mode=results$mode, depth=results$depth, sensitivity=format(results$sensitivity, digits=2), specificity = format(results$specificity, digits=2)))
+
+
+if ( onCMDLine ) dev.off()
+