From f131185e6857e947ef2db88b50756831628a813e Mon Sep 17 00:00:00 2001
From: Richard Leggett <richard.leggett@earlham.ac.uk>
Date: Fri, 3 Feb 2017 23:05:12 +0000
Subject: [PATCH] Error checking in R scripts

---
 Dockerfile                    |   4 +-
 bin/nanook_plot_lengths.R     |  60 ++--
 bin/nanook_plot_reference.R   | 550 ++++++++++++++++++----------------
 dist/NanoOK.jar               | Bin 338876 -> 338966 bytes
 src/nanook/NanoOK.java        |   2 +-
 src/nanook/NanoOKOptions.java |   3 +
 6 files changed, 335 insertions(+), 284 deletions(-)

diff --git a/Dockerfile b/Dockerfile
index 51fb4f2..19ca29c 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -19,8 +19,8 @@ RUN cd /usr ; git clone https://github.com/TGAC/NanoOK
 ENV NANOOK_DIR="/usr/NanoOK"
 RUN echo "export PATH=/usr/NanoOK/bin:${PATH}" >> ~/.bashrc
 RUN Rscript -e "install.packages('ggplot2', repos='https://cran.ma.imperial.ac.uk/')"
-RUN Rscript -e "install.packages('ggmap’, repos='https://cran.ma.imperial.ac.uk/')"
-RUN Rscript -e "install.packages('plyr’, repos='https://cran.ma.imperial.ac.uk/')"
+RUN Rscript -e "install.packages('ggmap', repos='https://cran.ma.imperial.ac.uk/')"
+RUN Rscript -e "install.packages('plyr', repos='https://cran.ma.imperial.ac.uk/')"
 RUN Rscript -e "install.packages('scales', repos='https://cran.ma.imperial.ac.uk/')"
 RUN Rscript -e "install.packages('gridExtra', repos='https://cran.ma.imperial.ac.uk/')"
 RUN Rscript -e "install.packages('reshape', repos='https://cran.ma.imperial.ac.uk/')"
\ No newline at end of file
diff --git a/bin/nanook_plot_lengths.R b/bin/nanook_plot_lengths.R
index f7a12fd..7c237c2 100755
--- a/bin/nanook_plot_lengths.R
+++ b/bin/nanook_plot_lengths.R
@@ -35,37 +35,53 @@ for (t in 1:3) {
     
     # Count vs length
     filename_lengths <- paste(analysisdir, "/", "all_",type,"_lengths.txt", sep="");
+    filename_kmers <- paste(analysisdir, "/", "all_",type,"_kmers.txt", sep="");
     
     if (file.exists(filename_lengths)) {
         data_lengths = read.table(filename_lengths, col.name=c("name", "length"))
 
-        if (format=="png") {
-            lengths_png <- paste(graphsdir, "/", "all_",type,"_lengths.png", sep="");
-            png(lengths_png, width=1200, height=800)
-            print(ggplot(data_lengths, aes(x=data_lengths$length), xlab="Length") + geom_histogram(binwidth=1000, fill=colourcode) + xlab("Length") +ylab("Count") + scale_x_continuous(limits=c(0, 35000)) + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+        if (nrow(data_lengths) > 1) {
+            if (format=="png") {
+                lengths_png <- paste(graphsdir, "/", "all_",type,"_lengths.png", sep="");
+                png(lengths_png, width=1200, height=800)
+                print(ggplot(data_lengths, aes(x=data_lengths$length), xlab="Length") + geom_histogram(binwidth=1000, fill=colourcode) + xlab("Length") +ylab("Count") + scale_x_continuous(limits=c(0, 35000)) + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            } else {
+                lengths_pdf <- paste(graphsdir, "/", "all_",type,"_lengths.pdf", sep="");
+                pdf(lengths_pdf, width=6, height=4)
+                print(ggplot(data_lengths, aes(x=data_lengths$length), xlab="Length") + geom_histogram(binwidth=1000, fill=colourcode) + xlab("Length") +ylab("Count") + scale_x_continuous(limits=c(0, 35000)) + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
         } else {
-            lengths_pdf <- paste(graphsdir, "/", "all_",type,"_lengths.pdf", sep="");
-            pdf(lengths_pdf, width=6, height=4)
-            print(ggplot(data_lengths, aes(x=data_lengths$length), xlab="Length") + geom_histogram(binwidth=1000, fill=colourcode) + xlab("Length") +ylab("Count") + scale_x_continuous(limits=c(0, 35000)) + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            cat("WARNING: No data in ", filename_lengths, "\n");
         }
-        garbage <- dev.off()
-
-        # Number of perfect 21mers verses length scatter
-        filename_kmers <- paste(analysisdir, "/", "all_",type,"_kmers.txt", sep="");
-        data_alignments = read.table(filename_kmers, header=TRUE)
+    } else {
+        cat("WARNING: Couldn't find", filename_lengths, "\n");
+    }
+    
+    # Number of perfect 21mers verses length scatter
+    if (file.exists(filename_kmers)) {
+        data_kmers = try(read.table(filename_kmers, header=TRUE), silent=TRUE)
 
-        if (format=="png") {
-            kmers_png <- paste(graphsdir, "/", "all_",type,"_21mers.png", sep="");
-            png(kmers_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$Length, y=data_alignments$nk21), xlab="Read length") + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Number of perfect 21mers") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(breaks=seq(0, 40000, 4000)) + scale_y_continuous(breaks=seq(0, 400, 20)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-            #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+        if (inherits(data_kmers, "try-error")) {
+            cat("WARNING: Couldn't read", filename_kmers,"\n");
         } else {
-            kmers_pdf <- paste(graphsdir, "/", "all_",type,"_21mers.pdf", sep="");
-            pdf(kmers_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$Length, y=data_alignments$nk21), xlab="Read length") + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Number of perfect 21mers") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(breaks=seq(0, 40000, 4000)) + scale_y_continuous(breaks=seq(0, 400, 20)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            if (nrow(data_kmers) > 1) {
+                if (format=="png") {
+                    kmers_png <- paste(graphsdir, "/", "all_",type,"_21mers.png", sep="");
+                    png(kmers_png, width=1200, height=800)
+                    print(ggplot(data_kmers, aes(x=data_kmers$Length, y=data_kmers$nk21), xlab="Read length") + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Number of perfect 21mers") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(breaks=seq(0, 40000, 4000)) + scale_y_continuous(breaks=seq(0, 400, 20)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+                    #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+                } else {
+                    kmers_pdf <- paste(graphsdir, "/", "all_",type,"_21mers.pdf", sep="");
+                    pdf(kmers_pdf, width=6, height=4)
+                    print(ggplot(data_kmers, aes(x=data_kmers$Length, y=data_kmers$nk21), xlab="Read length") + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Number of perfect 21mers") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(breaks=seq(0, 40000, 4000)) + scale_y_continuous(breaks=seq(0, 400, 20)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+                }
+                garbage <- dev.off()
+            } else {
+                cat("WARNING: No data in ", filename_kmers, "\n");
+            }
         }
-        garbage <- dev.off()
     } else {
-        cat("WARNING: Couldn't find", filename_lengths, "\n");
+        cat("WARNING: Couldn't find", filename_kmers, "\n");
     }
 }
diff --git a/bin/nanook_plot_reference.R b/bin/nanook_plot_reference.R
index 088546d..43ef583 100755
--- a/bin/nanook_plot_reference.R
+++ b/bin/nanook_plot_reference.R
@@ -38,19 +38,29 @@ if (format=="png") {
 
 # Plot GC% vs position
 data_gc_filename <- paste(analysisdir, "/", refid, "/", refid, "_gc.txt", sep="");
-data_gc = read.table(data_gc_filename, col.name=c("Position", "Coverage"))
-if (format=="png") {
-    png_gc <- paste(graphsdir, "/", refid, "/", refid, "_gc.png", sep="");
-    cat("Writing", png_gc, "\n");
-    png(png_gc, width=1600, height=400)
-    print(ggplot(data_gc, aes(x=data_gc$Position, y=data_gc$Coverage)) + geom_line(color="black") + ggtitle("GC content") + theme(text = element_text(size=textsize)) + xlab("Position") + ylab("GC %") + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+
+if (file.exists(data_gc_filename)) {
+    data_gc = read.table(data_gc_filename, col.name=c("Position", "Coverage"))
+    
+    if (nrow(data_gc) > 1) {
+        if (format=="png") {
+            png_gc <- paste(graphsdir, "/", refid, "/", refid, "_gc.png", sep="");
+            cat("Writing", png_gc, "\n");
+            png(png_gc, width=1600, height=400)
+            print(ggplot(data_gc, aes(x=data_gc$Position, y=data_gc$Coverage)) + geom_line(color="black") + ggtitle("GC content") + theme(text = element_text(size=textsize)) + xlab("Position") + ylab("GC %") + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+        } else {
+            pdf_gc <- paste(graphsdir, "/", refid, "/", refid, "_gc.pdf", sep="");
+            cat("Writing", pdf_gc, "\n");
+            pdf(pdf_gc, width=16, height=4)
+            print(ggplot(data_gc, aes(x=data_gc$Position, y=data_gc$Coverage)) + geom_line(color="black") + ggtitle("GC content") + theme(text = element_text(size=textsize)) + xlab("Position") + ylab("GC %") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+        }
+        garbage <- dev.off()
+    } else {
+        cat("WARNING: No data in ", data_gc_filename, "\n");
+    }
 } else {
-    pdf_gc <- paste(graphsdir, "/", refid, "/", refid, "_gc.pdf", sep="");
-    cat("Writing", pdf_gc, "\n");
-    pdf(pdf_gc, width=16, height=4)
-    print(ggplot(data_gc, aes(x=data_gc$Position, y=data_gc$Coverage)) + geom_line(color="black") + ggtitle("GC content") + theme(text = element_text(size=textsize)) + xlab("Position") + ylab("GC %") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+    cat("WARNING: Couldn't find", data_gc_filename, "\n");
 }
-garbage <- dev.off()
 
 listOfDataFrames <- NULL;
 count <-c(1);
@@ -77,6 +87,8 @@ for (t in 1:3) {
                     }
                 }
             }
+        } else {
+            cat("WARNING: No data in ", data_perfect_cumulative_filename, "\n");
         }
     }
 }
@@ -96,18 +108,22 @@ for (t in 1:3) {
 
     if (file.exists(data_coverage_filename)) {
         data_coverage = read.table(data_coverage_filename, col.name=c("Position", "Coverage"))
-        if (format=="png") {
-            png_coverage <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_coverage.png", sep="");
-            cat("Writing", png_coverage, "\n");
-            png(png_coverage, width=1600, height=400)
-            print(ggplot(data_coverage, aes(x=data_coverage$Position, y=data_coverage$Coverage)) + geom_line(color=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Position") + ylab("Mean coverage") + expand_limits(y = 0) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+        if (nrow(data_coverage) > 0) {
+            if (format=="png") {
+                png_coverage <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_coverage.png", sep="");
+                cat("Writing", png_coverage, "\n");
+                png(png_coverage, width=1600, height=400)
+                print(ggplot(data_coverage, aes(x=data_coverage$Position, y=data_coverage$Coverage)) + geom_line(color=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Position") + ylab("Mean coverage") + expand_limits(y = 0) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            } else {
+                pdf_coverage <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_coverage.pdf", sep="");
+                cat("Writing", pdf_coverage, "\n");
+                pdf(pdf_coverage, width=16, height=4)
+                print(ggplot(data_coverage, aes(x=data_coverage$Position, y=data_coverage$Coverage)) + geom_line(color=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Position") + ylab("Mean coverage") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + expand_limits(y = 0) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
         } else {
-            pdf_coverage <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_coverage.pdf", sep="");
-            cat("Writing", pdf_coverage, "\n");
-            pdf(pdf_coverage, width=16, height=4)
-            print(ggplot(data_coverage, aes(x=data_coverage$Position, y=data_coverage$Coverage)) + geom_line(color=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Position") + ylab("Mean coverage") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + expand_limits(y = 0) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            cat("WARNING: No data in ", data_coverage_filename, "\n");
         }
-        garbage <- dev.off()
     } else {
         cat("WARNING: Couldn't find", data_coverage_filename, "\n");
     }
@@ -118,19 +134,22 @@ for (t in 1:3) {
     if (file.exists(data_perfect_cumulative_filename)) {
         cat("Reading", data_perfect_cumulative_filename, "\n");
         data_perfect_cumulative = read.table(data_perfect_cumulative_filename, col.name=c("Size", "n", "Perfect"))
-
-        if (format=="png") {
-            png_perfect_cumulative <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_cumulative_perfect_kmers.png", sep="");
-            cat("Writing", png_perfect_cumulative, "\n");
-            png(png_perfect_cumulative, width=1200, height=800)
-            print(ggplot(data_perfect_cumulative, aes(x=data_perfect_cumulative$Size, y=data_perfect_cumulative$Perfect)) + geom_bar(stat="identity", width=0.7, fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("kmer size") + ylab("% reads with perfect kmer") + scale_x_continuous(limits=c(0, cum_maxk), breaks=seq(0,cum_maxk,by=50)) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+        if (nrow(data_perfect_cumulative) > 0) {
+            if (format=="png") {
+                png_perfect_cumulative <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_cumulative_perfect_kmers.png", sep="");
+                cat("Writing", png_perfect_cumulative, "\n");
+                png(png_perfect_cumulative, width=1200, height=800)
+                print(ggplot(data_perfect_cumulative, aes(x=data_perfect_cumulative$Size, y=data_perfect_cumulative$Perfect)) + geom_bar(stat="identity", width=0.7, fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("kmer size") + ylab("% reads with perfect kmer") + scale_x_continuous(limits=c(0, cum_maxk), breaks=seq(0,cum_maxk,by=50)) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            } else {
+                pdf_perfect_cumulative <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_cumulative_perfect_kmers.pdf", sep="");
+                cat("Writing", pdf_perfect_cumulative, "\n");
+                pdf(pdf_perfect_cumulative, width=6, height=4)
+                print(ggplot(data_perfect_cumulative, aes(x=data_perfect_cumulative$Size, y=data_perfect_cumulative$Perfect)) + geom_bar(stat="identity", width=0.7, fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("kmer size") + ylab("% reads with perfect kmer") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + scale_x_continuous(limits=c(0, cum_maxk), breaks=seq(0,cum_maxk,by=50)) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
         } else {
-            pdf_perfect_cumulative <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_cumulative_perfect_kmers.pdf", sep="");
-            cat("Writing", pdf_perfect_cumulative, "\n");
-            pdf(pdf_perfect_cumulative, width=6, height=4)
-            print(ggplot(data_perfect_cumulative, aes(x=data_perfect_cumulative$Size, y=data_perfect_cumulative$Perfect)) + geom_bar(stat="identity", width=0.7, fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("kmer size") + ylab("% reads with perfect kmer") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + scale_x_continuous(limits=c(0, cum_maxk), breaks=seq(0,cum_maxk,by=50)) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            cat("WARNING: No data in ", data_perfect_cumulative_filename, "\n");
         }
-        garbage <- dev.off()
     } else {
         cat("WARNING: Couldn't find", data_perfect_cumulative_filename, "\n");
     }
@@ -156,20 +175,23 @@ for (t in 1:3) {
     data_insertions_filename <- paste(analysisdir, "/", refid, "/", refid, "_",type,"_insertions.txt", sep="");
     cat("Reading", data_insertions_filename, "\n");
     if (file.exists(data_insertions_filename)) {
-        if (format=="png") {
-            png_insertions <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_insertions.png", sep="");
-            cat("Writing", png_insertions, "\n");
-            png(png_insertions, width=1200, height=800)
-            data_insertions = read.table(data_insertions_filename, col.name=c("Size", "Percent"))
-            print(ggplot(data_insertions, aes(x=data_insertions$Size, y=data_insertions$Percent)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Insertion size") + ylab("%") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+        data_insertions = read.table(data_insertions_filename, col.name=c("Size", "Percent"))
+        if (nrow(data_insertions) > 0) {
+            if (format=="png") {
+                png_insertions <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_insertions.png", sep="");
+                cat("Writing", png_insertions, "\n");
+                png(png_insertions, width=1200, height=800)
+                print(ggplot(data_insertions, aes(x=data_insertions$Size, y=data_insertions$Percent)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Insertion size") + ylab("%") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            } else {
+                pdf_insertions <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_insertions.pdf", sep="");
+                cat("Writing", pdf_insertions, "\n");
+                pdf(pdf_insertions, width=6, height=4)
+                print(ggplot(data_insertions, aes(x=data_insertions$Size, y=data_insertions$Percent)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Insertion size") + ylab("%") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
         } else {
-            pdf_insertions <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_insertions.pdf", sep="");
-            cat("Writing", pdf_insertions, "\n");
-            pdf(pdf_insertions, width=6, height=4)
-            data_insertions = read.table(data_insertions_filename, col.name=c("Size", "Percent"))
-            print(ggplot(data_insertions, aes(x=data_insertions$Size, y=data_insertions$Percent)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Insertion size") + ylab("%") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            cat("WARNING: No data in ", data_insertions_filename, "\n");
         }
-        garbage <- dev.off()
     } else {
         cat("WARNING: Couldn't find", data_insertions_filename, "\n");
     }
@@ -178,20 +200,23 @@ for (t in 1:3) {
     data_deletions_filename <- paste(analysisdir, "/", refid, "/", refid, "_",type,"_deletions.txt", sep="");
     cat("Reading", data_deletions_filename, "\n");
     if (file.exists(data_deletions_filename)) {
-        if (format=="png") {
-            png_deletions <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_deletions.png", sep="");
-            cat("Writing", png_deletions, "\n");
-            png(png_deletions, width=1200, height=800)
-            data_deletions = read.table(data_deletions_filename, col.name=c("Size", "Percent"))
-            print(ggplot(data_deletions, aes(x=data_deletions$Size, y=data_deletions$Percent)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Deletion size") + ylab("%") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+        data_deletions = read.table(data_deletions_filename, col.name=c("Size", "Percent"))
+        if (nrow(data_deletions) > 0) {
+            if (format=="png") {
+                png_deletions <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_deletions.png", sep="");
+                cat("Writing", png_deletions, "\n");
+                png(png_deletions, width=1200, height=800)
+                print(ggplot(data_deletions, aes(x=data_deletions$Size, y=data_deletions$Percent)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Deletion size") + ylab("%") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            } else {
+                pdf_deletions <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_deletions.pdf", sep="");
+                cat("Writing", pdf_deletions, "\n");
+                pdf(pdf_deletions, width=6, height=4)
+                print(ggplot(data_deletions, aes(x=data_deletions$Size, y=data_deletions$Percent)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Deletion size") + ylab("%") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
         } else {
-            pdf_deletions <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_deletions.pdf", sep="");
-            cat("Writing", pdf_deletions, "\n");
-            pdf(pdf_deletions, width=6, height=4)
-            data_deletions = read.table(data_deletions_filename, col.name=c("Size", "Percent"))
-            print(ggplot(data_deletions, aes(x=data_deletions$Size, y=data_deletions$Percent)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Deletion size") + ylab("%") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            cat("WARNING: No data in ", data_deletions_filename, "\n");
         }
-        garbage <- dev.off()
     } else {
         cat("WARNING: Couldn't find", data_deletions_filename, "\n");
     }
@@ -203,200 +228,203 @@ for (t in 1:3) {
     
     if (file.exists(input_filename)) {
         data_alignments = read.table(input_filename, header=TRUE);
-
-        # Length vs Identity histograms
-        if (format=="png") {
-            identity_hist_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_hist.png", sep="")
-            cat("Writing", identity_hist_png, "\n");
-            png(identity_hist_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryPercentIdentity)) + geom_histogram(fill=colourcode) + xlab("Read identity %") +ylab("Count") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) )
-        } else {
-            identity_hist_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_hist.pdf", sep="")
-            cat("Writing", identity_hist_pdf, "\n");
-            pdf(identity_hist_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryPercentIdentity)) + geom_histogram(fill=colourcode) + xlab("Read identity %") +ylab("Count") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-        }
-        garbage <- dev.off()
-        
-        # GC histogram
-        if (format=="png") {
-            identity_hist_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_GC_hist.png", sep="")
-            cat("Writing", identity_hist_png, "\n");
-            png(identity_hist_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryGC)) + geom_histogram(fill=colourcode, binwidth=1) + xlab("GC %") +ylab("Read count") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) + scale_x_continuous(limits=c(0, 100)) )
-        } else {
-            identity_hist_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_GC_hist.pdf", sep="")
-            cat("Writing", identity_hist_pdf, "\n");
-            pdf(identity_hist_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryGC)) + geom_histogram(fill=colourcode, binwidth = 1) + xlab("GC %") +ylab("Read count") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) + scale_x_continuous(limits=c(0, 100)))
-        }
-        garbage <- dev.off()
-
-        # Identity vs Length Scatter plots
-        if (format=="png") {
-            identity_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_scatter.png", sep="");
-            cat("Writing", identity_scatter_pdf, "\n");
-            png(identity_scatter_pdf, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-            #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
-        } else {
-            identity_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_scatter.pdf", sep="");
-            cat("Writing", identity_scatter_pdf, "\n");
-            pdf(identity_scatter_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-        }
-        garbage <- dev.off()
-
-        # Identity vs Length heatmap
-        if (format=="png") {
-            identity_heatmap_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_heatmap.png", sep="");
-            cat("Writing", identity_heatmap_png, "\n");
-            png(identity_heatmap_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_bin2d(drop=TRUE, binwidth=c(500,2)) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=(textsize*0.75))) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) + scale_fill_gradientn(colours=rev(rainbow(n=30, end=4/6)), na.value='blue'))
-            #limits=c(0,50), breaks=seq(0, 40, by=10), colours=rainbow(4)
-        } else {
-            identity_heatmap_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_heatmap.pdf", sep="");
-            cat("Writing", identity_heatmap_pdf, "\n");
-            pdf(identity_heatmap_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_bin2d(drop=TRUE, binwidth=c(500,2)) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=(textsize*0.75))) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust))+ scale_fill_gradientn(colours=rev(rainbow(n=30, end=4/6)), na.value='blue'))
-        }
-        garbage <- dev.off()
-
-        # Identity vs Length heatmap zoomed
-        if (format=="png") {
-            identity_heatmap_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_heatmap_zoom.png", sep="");
-            cat("Writing", identity_heatmap_png, "\n");
-            png(identity_heatmap_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_bin2d(drop=TRUE, binwidth=c(500,1)) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=(textsize*0.75))) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) + scale_fill_gradientn(colours=rev(rainbow(n=30, end=4/6)), na.value='blue')+ scale_y_continuous(limits=c(60, 100)))
-            #limits=c(0,50), breaks=seq(0, 40, by=10), colours=rainbow(4)
-        } else {
-            identity_heatmap_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_heatmap_zoom.pdf", sep="");
-            cat("Writing", identity_heatmap_pdf, "\n");
-            pdf(identity_heatmap_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_bin2d(drop=TRUE, binwidth=c(500,1)) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=(textsize*0.75))) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust))+ scale_fill_gradientn(colours=rev(rainbow(n=30, end=4/6)), na.value='blue')+ scale_y_continuous(limits=c(60, 100)))
-        }
-        garbage <- dev.off()
-
-        # Identity boxplot
-        listOfDataFrames[[count]] <- data.frame(Readset=type, Variable=data_alignments$QueryPercentIdentity);
-        count <- count + 1;
-
-        # Alignment identity vs. Fraction of read aligned scatter plots
-        if (format=="png") {
-            aid_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_read_fraction_vs_alignment_identity_scatter.png", sep="");
-            cat("Writing", aid_scatter_png, "\n");
-            png(aid_scatter_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$PercentQueryAligned, y=data_alignments$AlignmentPercentIdentity)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Percentage of read aligned") +ylab("Alignment identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 105)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) )
-            #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
-        } else {
-            aid_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_read_fraction_vs_alignment_identity_scatter.pdf", sep="");
-            cat("Writing", aid_scatter_pdf, "\n");
-            pdf(aid_scatter_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$PercentQueryAligned, y=data_alignments$AlignmentPercentIdentity)) + geom_point(shape=pointshape, alpha=0.4, color=colourcode) + xlab("Percentage of read aligned") +ylab("Alignment identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 105)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-        }
-        garbage <- dev.off()
-
-        # Query identity vs. Fraction of read aligned scatter plots
-        if (format=="png") {
-            qid_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_read_fraction_vs_query_identity_scatter.png", sep="");
-            cat("Writing", qid_scatter_png, "\n");
-            png(qid_scatter_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$PercentQueryAligned, y=data_alignments$QueryPercentIdentity)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Percentage of read aligned") +ylab("Alignment identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 105)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-            #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
-        } else {
-            qid_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_read_fraction_vs_query_identity_scatter.pdf", sep="");
-            cat("Writing", qid_scatter_pdf, "\n");
-            pdf(qid_scatter_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$PercentQueryAligned, y=data_alignments$QueryPercentIdentity)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Percentage of read aligned") +ylab("Alignment identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 105)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-        }
-        garbage <- dev.off()
-
-        # Best perfect sequence vs. length scatters
-        if (format=="png") {
-            best_perf_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_longest_perfect_vs_length_scatter.png", sep="");
-            cat("Writing", best_perf_scatter_png, "\n");
-            png(best_perf_scatter_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$LongestPerfectKmer)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Longest perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-            #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
-        } else {
-            best_perf_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_longest_perfect_vs_length_scatter.pdf", sep="");
-            cat("Writing", best_perf_scatter_pdf, "\n");
-            pdf(best_perf_scatter_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$LongestPerfectKmer)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Longest perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-        }
-        garbage <- dev.off()
-
-        # Best perfect sequence vs. length scatters zoomed
-        if (format=="png") {
-            best_perf_zoom_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_longest_perfect_vs_length_zoom_scatter.png", sep="");
-            cat("Writing", best_perf_zoom_scatter_png, "\n");
-            png(best_perf_zoom_scatter_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$LongestPerfectKmer)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Longest perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 10000)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-            #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
-        } else {
-            best_perf_zoom_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_longest_perfect_vs_length_zoom_scatter.pdf", sep="");
-            cat("Writing", best_perf_zoom_scatter_pdf, "\n");
-            pdf(best_perf_zoom_scatter_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$LongestPerfectKmer)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Longest perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 10000)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-        }
-        garbage <- dev.off()
-
-        # Plot %reads vs best perfect kmer
-        cat(data_alignments$LongestPerfectKmer);
-        
-        hdf <- hist(breaks=seq(0,maxk,by=10), x=data_alignments$LongestPerfectKmer, plot=FALSE, right=FALSE); # bins are 0-9, 10-19, 20-29 etc.
-        hdf$density = hdf$counts/sum(hdf$counts)*100
-        tdf <- data.frame(Pos=hdf$mids, Counts=hdf$density);
-
-        if (format=="png") {
-            png_perfect_best <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_best_perfect_kmers.png", sep="");
-            cat("Writing", png_perfect_best, "\n");
-            png(png_perfect_best, width=1200, height=800)
-            print(ggplot(tdf, aes(Pos, Counts)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Best perfect kmer") + ylab("% reads") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + scale_x_continuous(limits=c(0, cum_maxk), breaks=seq(0,cum_maxk,by=50)) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-        } else {
-            pdf_perfect_best <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_best_perfect_kmers.pdf", sep="");
-            cat("Writing", pdf_perfect_best, "\n");
-            pdf(pdf_perfect_best, width=6, height=4)
-            print(ggplot(tdf, aes(Pos, Counts)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Best perfect kmer") + ylab("% reads") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + scale_x_continuous(limits=c(0, cum_maxk), breaks=seq(0,cum_maxk,by=50)) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-        }
-        garbage <- dev.off()
-
-
-        # Number of perfect 21mers verses length scatter
-        if (format=="png") {
-            nk21_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_nk21_vs_length_scatter.png", sep="");
-            cat("Writing", nk21_scatter_png, "\n");
-            png(nk21_scatter_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$nk21)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Number of perfect 21mers") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-            #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
-        } else {
-            nk21_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_nk21_vs_length_scatter.pdf", sep="");
-            cat("Writing", nk21_scatter_pdf, "\n");
-            pdf(nk21_scatter_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$nk21)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Number of perfect 21mers") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-        }
-        garbage <- dev.off()
-
-        # Mean perfect sequence vs. length scatters
-        #mean_perf_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_mean_perfect_vs_length_scatter.pdf", sep="");
-        #pdf(mean_perf_scatter_pdf, height=4, width=6)
-        #print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$MeanPerfectKmer), xlab="Read length") + geom_point(shape=pointshape, alpha=pointalpha) + xlab("Read length") +ylab("Mean perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 10000)))
-        #garbage <- dev.off()
-
-        # Percentage of read aligned vs read length
-        if (format=="png") {
-            output_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_percent_aligned_vs_length_scatter.png", sep="");
-            cat("Writing", output_png, "\n");
-            png(output_png, width=1200, height=800)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$PercentQueryAligned)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Percentage of read aligned") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
-            #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+        if (nrow(data_alignments) > 1) {
+            # Length vs Identity histograms
+            if (format=="png") {
+                identity_hist_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_hist.png", sep="")
+                cat("Writing", identity_hist_png, "\n");
+                png(identity_hist_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryPercentIdentity)) + geom_histogram(fill=colourcode) + xlab("Read identity %") +ylab("Count") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) )
+            } else {
+                identity_hist_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_hist.pdf", sep="")
+                cat("Writing", identity_hist_pdf, "\n");
+                pdf(identity_hist_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryPercentIdentity)) + geom_histogram(fill=colourcode) + xlab("Read identity %") +ylab("Count") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
+            
+            # GC histogram
+            if (format=="png") {
+                identity_hist_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_GC_hist.png", sep="")
+                cat("Writing", identity_hist_png, "\n");
+                png(identity_hist_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryGC)) + geom_histogram(fill=colourcode, binwidth=1) + xlab("GC %") +ylab("Read count") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) + scale_x_continuous(limits=c(0, 100)) )
+            } else {
+                identity_hist_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_GC_hist.pdf", sep="")
+                cat("Writing", identity_hist_pdf, "\n");
+                pdf(identity_hist_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryGC)) + geom_histogram(fill=colourcode, binwidth = 1) + xlab("GC %") +ylab("Read count") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) + scale_x_continuous(limits=c(0, 100)))
+            }
+            garbage <- dev.off()
+
+            # Identity vs Length Scatter plots
+            if (format=="png") {
+                identity_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_scatter.png", sep="");
+                cat("Writing", identity_scatter_pdf, "\n");
+                png(identity_scatter_pdf, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+                #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+            } else {
+                identity_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_scatter.pdf", sep="");
+                cat("Writing", identity_scatter_pdf, "\n");
+                pdf(identity_scatter_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
+
+            # Identity vs Length heatmap
+            if (format=="png") {
+                identity_heatmap_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_heatmap.png", sep="");
+                cat("Writing", identity_heatmap_png, "\n");
+                png(identity_heatmap_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_bin2d(drop=TRUE, binwidth=c(500,2)) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=(textsize*0.75))) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) + scale_fill_gradientn(colours=rev(rainbow(n=30, end=4/6)), na.value='blue'))
+                #limits=c(0,50), breaks=seq(0, 40, by=10), colours=rainbow(4)
+            } else {
+                identity_heatmap_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_heatmap.pdf", sep="");
+                cat("Writing", identity_heatmap_pdf, "\n");
+                pdf(identity_heatmap_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_bin2d(drop=TRUE, binwidth=c(500,2)) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=(textsize*0.75))) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust))+ scale_fill_gradientn(colours=rev(rainbow(n=30, end=4/6)), na.value='blue'))
+            }
+            garbage <- dev.off()
+
+            # Identity vs Length heatmap zoomed
+            if (format=="png") {
+                identity_heatmap_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_heatmap_zoom.png", sep="");
+                cat("Writing", identity_heatmap_png, "\n");
+                png(identity_heatmap_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_bin2d(drop=TRUE, binwidth=c(500,1)) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=(textsize*0.75))) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) + scale_fill_gradientn(colours=rev(rainbow(n=30, end=4/6)), na.value='blue')+ scale_y_continuous(limits=c(60, 100)))
+                #limits=c(0,50), breaks=seq(0, 40, by=10), colours=rainbow(4)
+            } else {
+                identity_heatmap_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_length_vs_identity_heatmap_zoom.pdf", sep="");
+                cat("Writing", identity_heatmap_pdf, "\n");
+                pdf(identity_heatmap_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$QueryPercentIdentity)) + geom_bin2d(drop=TRUE, binwidth=c(500,1)) + xlab("Length") +ylab("Read identity %") + ggtitle(type) + theme(text = element_text(size=(textsize*0.75))) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust))+ scale_fill_gradientn(colours=rev(rainbow(n=30, end=4/6)), na.value='blue')+ scale_y_continuous(limits=c(60, 100)))
+            }
+            garbage <- dev.off()
+
+            # Identity boxplot
+            listOfDataFrames[[count]] <- data.frame(Readset=type, Variable=data_alignments$QueryPercentIdentity);
+            count <- count + 1;
+
+            # Alignment identity vs. Fraction of read aligned scatter plots
+            if (format=="png") {
+                aid_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_read_fraction_vs_alignment_identity_scatter.png", sep="");
+                cat("Writing", aid_scatter_png, "\n");
+                png(aid_scatter_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$PercentQueryAligned, y=data_alignments$AlignmentPercentIdentity)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Percentage of read aligned") +ylab("Alignment identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 105)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) )
+                #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+            } else {
+                aid_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_read_fraction_vs_alignment_identity_scatter.pdf", sep="");
+                cat("Writing", aid_scatter_pdf, "\n");
+                pdf(aid_scatter_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$PercentQueryAligned, y=data_alignments$AlignmentPercentIdentity)) + geom_point(shape=pointshape, alpha=0.4, color=colourcode) + xlab("Percentage of read aligned") +ylab("Alignment identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 105)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
+
+            # Query identity vs. Fraction of read aligned scatter plots
+            if (format=="png") {
+                qid_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_read_fraction_vs_query_identity_scatter.png", sep="");
+                cat("Writing", qid_scatter_png, "\n");
+                png(qid_scatter_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$PercentQueryAligned, y=data_alignments$QueryPercentIdentity)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Percentage of read aligned") +ylab("Alignment identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 105)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+                #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+            } else {
+                qid_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_read_fraction_vs_query_identity_scatter.pdf", sep="");
+                cat("Writing", qid_scatter_pdf, "\n");
+                pdf(qid_scatter_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$PercentQueryAligned, y=data_alignments$QueryPercentIdentity)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Percentage of read aligned") +ylab("Alignment identity %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 105)) + scale_y_continuous(limits=c(0, 100)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
+
+            # Best perfect sequence vs. length scatters
+            if (format=="png") {
+                best_perf_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_longest_perfect_vs_length_scatter.png", sep="");
+                cat("Writing", best_perf_scatter_png, "\n");
+                png(best_perf_scatter_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$LongestPerfectKmer)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Longest perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+                #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+            } else {
+                best_perf_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_longest_perfect_vs_length_scatter.pdf", sep="");
+                cat("Writing", best_perf_scatter_pdf, "\n");
+                pdf(best_perf_scatter_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$LongestPerfectKmer)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Longest perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
+
+            # Best perfect sequence vs. length scatters zoomed
+            if (format=="png") {
+                best_perf_zoom_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_longest_perfect_vs_length_zoom_scatter.png", sep="");
+                cat("Writing", best_perf_zoom_scatter_png, "\n");
+                png(best_perf_zoom_scatter_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$LongestPerfectKmer)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Longest perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 10000)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+                #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+            } else {
+                best_perf_zoom_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_longest_perfect_vs_length_zoom_scatter.pdf", sep="");
+                cat("Writing", best_perf_zoom_scatter_pdf, "\n");
+                pdf(best_perf_zoom_scatter_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$LongestPerfectKmer)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Longest perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 10000)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
+
+            # Plot %reads vs best perfect kmer
+            cat(data_alignments$LongestPerfectKmer);
+            
+            hdf <- hist(breaks=seq(0,maxk,by=10), x=data_alignments$LongestPerfectKmer, plot=FALSE, right=FALSE); # bins are 0-9, 10-19, 20-29 etc.
+            hdf$density = hdf$counts/sum(hdf$counts)*100
+            tdf <- data.frame(Pos=hdf$mids, Counts=hdf$density);
+
+            if (format=="png") {
+                png_perfect_best <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_best_perfect_kmers.png", sep="");
+                cat("Writing", png_perfect_best, "\n");
+                png(png_perfect_best, width=1200, height=800)
+                print(ggplot(tdf, aes(Pos, Counts)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Best perfect kmer") + ylab("% reads") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + scale_x_continuous(limits=c(0, cum_maxk), breaks=seq(0,cum_maxk,by=50)) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            } else {
+                pdf_perfect_best <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_best_perfect_kmers.pdf", sep="");
+                cat("Writing", pdf_perfect_best, "\n");
+                pdf(pdf_perfect_best, width=6, height=4)
+                print(ggplot(tdf, aes(Pos, Counts)) + geom_bar(stat="identity", fill=colourcode) + ggtitle(type) + theme(text = element_text(size=textsize)) + xlab("Best perfect kmer") + ylab("% reads") + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.y=element_text(vjust=0.2)) + theme(axis.title.x=element_text(vjust=-0.2)) + scale_x_continuous(limits=c(0, cum_maxk), breaks=seq(0,cum_maxk,by=50)) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
+
+
+            # Number of perfect 21mers verses length scatter
+            if (format=="png") {
+                nk21_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_nk21_vs_length_scatter.png", sep="");
+                cat("Writing", nk21_scatter_png, "\n");
+                png(nk21_scatter_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$nk21)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Number of perfect 21mers") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+                #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+            } else {
+                nk21_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_nk21_vs_length_scatter.pdf", sep="");
+                cat("Writing", nk21_scatter_pdf, "\n");
+                pdf(nk21_scatter_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$nk21)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Number of perfect 21mers") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
+
+            # Mean perfect sequence vs. length scatters
+            #mean_perf_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_mean_perfect_vs_length_scatter.pdf", sep="");
+            #pdf(mean_perf_scatter_pdf, height=4, width=6)
+            #print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$MeanPerfectKmer), xlab="Read length") + geom_point(shape=pointshape, alpha=pointalpha) + xlab("Read length") +ylab("Mean perfect kmer") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 10000)))
+            #garbage <- dev.off()
+
+            # Percentage of read aligned vs read length
+            if (format=="png") {
+                output_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_percent_aligned_vs_length_scatter.png", sep="");
+                cat("Writing", output_png, "\n");
+                png(output_png, width=1200, height=800)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$PercentQueryAligned)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Percentage of read aligned") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+                #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+            } else {
+                output_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_percent_aligned_vs_length_scatter.pdf", sep="");
+                cat("Writing", output_pdf, "\n");
+                pdf(output_pdf, width=6, height=4)
+                print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$PercentQueryAligned)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Percentage of read aligned") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            }
+            garbage <- dev.off()
         } else {
-            output_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_percent_aligned_vs_length_scatter.pdf", sep="");
-            cat("Writing", output_pdf, "\n");
-            pdf(output_pdf, width=6, height=4)
-            print(ggplot(data_alignments, aes(x=data_alignments$QueryLength, y=data_alignments$PercentQueryAligned)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Read length") +ylab("Percentage of read aligned") + ggtitle(type) + theme(text = element_text(size=textsize)) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)))
+            cat("WARNING: No data in ", input_filename, "\n");
         }
-        garbage <- dev.off()
     } else {
         cat("WARNING: Couldn't find", input_filename, "\n");
     }
@@ -409,19 +437,23 @@ for (t in 1:3) {
 
     if (file.exists(input_kmers)) {
         data_kmers = read.table(input_kmers, header=TRUE);
-        if (format=="png") {
-            kmer_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_kmer_scatter.png", sep="");
-            cat("Writing", kmer_scatter_png, "\n");
-            png(kmer_scatter_png, width=1200, height=1200)
-            print(ggplot(data_kmers, aes(x=data_kmers$RefPc, y=data_kmers$ReadPc)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Reference abundance %") +ylab("Reads abundance %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 0.3)) + scale_y_continuous(limits=c(0, 0.3)) + geom_text(aes(label=data_kmers$Kmer), size=4) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) )
-            #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+        if (nrow(data_kmers) > 1) {
+            if (format=="png") {
+                kmer_scatter_png <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_kmer_scatter.png", sep="");
+                cat("Writing", kmer_scatter_png, "\n");
+                png(kmer_scatter_png, width=1200, height=1200)
+                print(ggplot(data_kmers, aes(x=data_kmers$RefPc, y=data_kmers$ReadPc)) + geom_point(shape=pointshape, size=pointsize, alpha=pointalpha, color=colourcode) + xlab("Reference abundance %") +ylab("Reads abundance %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 0.3)) + scale_y_continuous(limits=c(0, 0.3)) + geom_text(aes(label=data_kmers$Kmer), size=4) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) )
+                #grid.edit("geom_point.points", grep = TRUE, gp = gpar(lwd = pointwidth))
+            } else {
+                kmer_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_kmer_scatter.pdf", sep="");
+                cat("Writing", kmer_scatter_pdf, "\n");
+                pdf(kmer_scatter_pdf, width=6, height=6)
+                print(ggplot(data_kmers, aes(x=data_kmers$RefPc, y=data_kmers$ReadPc)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Reference abundance %") +ylab("Reads abundance %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 0.3)) + scale_y_continuous(limits=c(0, 0.3)) + geom_text(aes(label=data_kmers$Kmer), size=1) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) )
+            }
+            garbage <- dev.off()
         } else {
-            kmer_scatter_pdf <- paste(graphsdir, "/", refid, "/", refid, "_",type,"_kmer_scatter.pdf", sep="");
-            cat("Writing", kmer_scatter_pdf, "\n");
-            pdf(kmer_scatter_pdf, width=6, height=6)
-            print(ggplot(data_kmers, aes(x=data_kmers$RefPc, y=data_kmers$ReadPc)) + geom_point(shape=pointshape, alpha=pointalpha, color=colourcode) + xlab("Reference abundance %") +ylab("Reads abundance %") + ggtitle(type) + theme(text = element_text(size=textsize)) + scale_x_continuous(limits=c(0, 0.3)) + scale_y_continuous(limits=c(0, 0.3)) + geom_text(aes(label=data_kmers$Kmer), size=1) + theme(plot.margin = unit(c(0.02,0.02,0.04,0.02), "npc")) + theme(axis.title.x=element_text(vjust=-xvjust)) + theme(axis.title.y=element_text(vjust=yvjust)) )
+            cat("WARNING: No data in ", input_kmers, "\n");
         }
-        garbage <- dev.off()
     } else {
         cat("WARNING: Couldn't find", input_kmers, "\n");
     }
diff --git a/dist/NanoOK.jar b/dist/NanoOK.jar
index 8af76d5bfc8b4cfb650a419fa0867415a397acfd..15dea865e3f5b19b690a8037ce2081433a2a10e9 100644
GIT binary patch
delta 11025
zcmb_CcYIaF(z~;JayB<7g^-Y5C?ODn2|)#U7zj;@^k)GP2#`=y2oXetD<~*kh_VVQ
z@+>@5R8$B)`?EjGvpx|~Ecie~ij87`Z)WdJ?#cE1zJI=tADP)bXUgvE?9A-m^ZaKC
z8}=kbj~F7Pk1@sACreA#MRyP%FD+TzeP?%GxL;v&@X+~?U0=PQuS$e2h5FpDd>`R1
z^JH@UufOr{+9nXuk%9A%^sQx#{!$rZy7iMhL7z{7*y0Oq*y!W4V3Q5^_?UouBUwJ&
zM_@C7`w47O@HfH_D0q<YLkb=yyjAHdQ{fgrwO=LhnhmcLeuKc91l}U>whjOCu`aNQ
zI`1fWm+*T8-dFH%!XFU$kibU-J|^%9flmqi$A-^*tXE)cd+5pM1ojg6B9i&@VV@0O
z`q+Q~#aHxvKY_0ad}DLd`PPQ-d~8@?6bI<d_cr|3$3_LX4k|e0gCF3q4L|zWRj`S`
zkx2F-{G{M0b&f^i>*ESe_~2*wMZvEMepB$f(x0Tkm4lT?P$Cf}h*nA@`FOlYRw6}-
zR3*}sNLM05iA*K3l*m@1HN9`+6KzF1C32K#uSBjA9hB(kWBIx^9oj}Jaj_EJ&}{nt
zbg1+zajg>LiEV-s6O|}cVv-V*kz7nsqRbZMKK7xQN?@82(+ST|VkY5P1S*KKQi&>`
zs1~yc)Zm=N99z^{=n?eK7b&h2b8Rut$G)ZMe<iL*D~b6^+@Snoff6?=u~3OcO5kX2
zQev?ZOGxU?O58$tsS?Ww-%9#!Q)0OiE0kEN#41~?w#6EHx7Nqb&?rxf+X<|*#d;;`
zeC)iqLy0?;*g!;DiMy1zTM6{cc(F-|dr<$q^yWSSn-yQCIBL6}dRqwm%@z+R@gR{;
zRpKEf9#&$j5`R}>n-Y)E<42WvOo{DEJZ_68d^}4$X^R~`-kR`JKHiq_(>|Ufo>bzQ
zD6v!QQlego-L`nv$2*CU{~%L5N8ot^FA#W9iI<dkS&3I@4F4qXDuLGsyiVW^CEg_b
z7J;`3;6}(7@1TYB*O`#fKVN)A4?af66rT|IREhu4c+uhW#U3R-x5ZvteBt9m(QN#3
zv5&x)MDdjp`<3__o4Pm)QZG^Buo6EiaYTuqNYYUyjwx|mi4#ivti&(2_*IGDe0-|7
z!xq2$_zZE<7N=}+T9<q(it6jL;IFoP9Sx~tbKx6aUVkDNVt}vTl?Q_aKT|)rGgJt^
zNbfF$arRyt7z(UU>jGK4@E||LS&Pcj%BrfF1*6K!%WBFhCzs7}_(Ac9!!^Il;dk?m
z4&OxJ9==h}D1vVOd-+Dr<gC&;wGO|J+RWjb2}jJCKc}{=!r}K5{T2d$Bk%x$2MIhx
z;9&w=IdcD<Z*%w~1Rf>u7=i5s9w+bwkvvJQ9R!{t@HF3+@9<}+xs$*y0`&xT6L^-u
zKL|WW;CWtmqB1$b;V+P=7YV#X;ANajDHkn!JNy;?PtHX3lyZl^O5=Ks*K%f+<K!-K
z`0MoK4dQr{z*{8qZ36$I?mK*&{>HKrQXKv+f6w9X6ZkiQ4+wn7->ZB1rbK-E5&ziX
zpAh(zz<&sQ#rIRRZ{>6N*ZdoYAK<lmtq{qn4*#D2mowk6l3~M#51G(!&?s9pIO420
zSGOZQEk&G{%#obTEhbH`bR>Xuq>xfSb|H+2upAiyDfMs>NITx<r#Q3AuCJ{rom}fk
zB_nlmPgoS?lTo&G92qTR^p>8`T`c2=^uC@jTf`#oh)ZD~7jZI57xaP$Ba+BRvNi7o
zsrr{*kfaVmiX&S<k&f>Tc_JPgy?R4I*OoMn1O%|wc-64$DkhcHWLK4EPp+ydo6{rP
zmWhsRC6gSP%+EMDj}*P7H{?g8Ix<bB>)1<0ek4CAGaSB;f2maoT$|e+BxjLVRF=+~
zKDDxJiUT?HAflqGc6#|72MS1KCJta#2@HyEjXmVga%46?RJWxfS)VC^e9;=`*0B$C
zEa^tPW@Sr*cs6Iez>#fa6gC#&-vtij(hS-gH4~~!Yvz<W(19Lx1gY13$`ke9`#?X@
zmY>nb`aqF5DBJ0Fec|eeGmgxG6uqP`j7mBv+dKSozQ^I85!g%M3%*CM{subggnrP<
zKM!pPy`Z-PQQ$aGNE&LOh~`#SGnHoA1&SPqrI*(uEq@!kGFNtRU=F?LObmj+#ni1O
zZ~=9%Bf5@gUzsO6;f(SfnJ+s#(2i!*$~UU4bjtL~sd`&KD9%18FL0!>VIe<*v^~&S
z`ZqdD7xWFgKa6b=F}k#(dR7^(Lw2@l>TCN$jwnKBTG}7hq6NCvmAsv(Gpl%nKDHKG
z)b||#i=}@S<Gg}h!%}%AufldUW-8#bSqh)SYn%GlVSlcxKM(uYhxO;X`ZwVF1+M;$
z*k9=CFT%c``u?E8oAAM67r_$j-|Xt&g8ikgJ}Fq{>XU+7UH#jzzueU)`W3D|(XYh5
zzj4N^@WE;q0SQ>+>XU%Au09F4-PI@hb*?_quXpu{zRuOZ1N(RKT7Q#?H=x-ulpFJ7
zG&dLny1^LJ4aSgeFotu3F`OHWf!tt><pyIoHy9(j!5GvH#<*@UhIWIuazB==s2C%>
z!5HKX#yD><hI)fB+8c}k-(ZaS24mPa7$d*I82k;!_-`;~0E00J7>s$qU`z!D@Ad~9
zm=VkaObP~LUNHFipd!o+rjN<NV9XB&V~Q{s^Mk>d91O<%U@)czgE2c8j0wVE%n=6r
zF;AEVrV4{GTOeLSsls5)76xO&Fc@=&!I(A-#-w2|<_&`}br_7<!(a?824kEt7^95A
z7-$IhpeSk{V5Tt`(+v8=x7Y^$9sXg)G0!pGATGmqSPmfkp?=;l3`_iQ{mmn=NbvuN
zZ(A{741>fG{uAPDY$QL5^XL2+zKdgU9>)weKW_6AjK$*fpZPB&3Evhlj4X()yl0s8
zBsNem!U*$l2`AxBfFK|o*?#4}1(}OctsgbBh`dfanYE*d-AChL%b3M(H6qOTSVT|+
z((~WZ<54p|$xpe21!Se|WS4}?T7j}wddV7YgkeNMF{k;NpfLBeE_tkJT!Cq=@toEl
zyn&eOkBx#nh1N~I1~!W>WPm&-?fcoq{477m1fCgF8s~wn^EhPAKqUdk^#XoOA07{*
z{1AVJ3{5$}V?<7vp|KY8Ff>@6%rW825B-_>N>m-;ZQ0IYmPK!|u*PfhBKInx&ltmk
z3XDXzm<X+MoBRXoKCfX%28o9u@stqpv!290&m~I=TQV#I|EwvJk~3LXzcms1@xwa0
z6tb=$x(m@+vhdIyO`Cw_G_^W4wHWIPF=m5xVGy7*o5f*Px@B1nv|UcHNXqa=!{SmX
zMh!2QLfT-`fG**wfw9g(xt*HHjbakzHi?Z34ZNcsHVM*D_KZpY&!gwiA`me2SCimn
z)UbInTxP;WG45~1o(Q_J!ybTmbZyQGaW}Wb7O<Xd4kQ?NFU6yMBDN%}WF^oFTb8j@
z7PHQ?8(0K|6#(~t#-gI4qtOGgGz~aEFGqIT+XtT0WmDh~j%8RGT-J>&83+vrn#QsR
zEstYqT#W{njIs*~3VNWxWJocq@opJZql`<-q3w{y)nIHw<F;fIAl0mVTh`Oe4A=mc
z5|~d=E}Uhf0%MLsxoJT=Jz5R}Q0~v=klT~wX5)5m$t0vh23jYLwSY{>!WOP_K$nqb
zGZ<?_b_*1QHLe-pqf?=co;MYuLz@Lgxm-9iI8+X8LL3Dy;<28@ZG*&ontL0<Bp0#O
zBpFLeGnV4!#O0pEShySF1=C<`ciOWt6zQOnDX3yOFx$m6$TfArcn(e*>v6*!JssSW
znCs#S%@`~AM&dRXM(ko^cCm3`Aq{Tu<itYXa6UU7#$6H`q@U%&O?d6(&Q?QQgGATM
zMMgRz81`W^pu?5K9*<thU^UKh4U2=@St8UKy|nDn7%K9N6ET()v}q!4Y4gl1dnALn
z2Fwf?CpdHp%6@kSw0FtA8)a`~v2YJ-37eZ_J19FpC_6DIyCupt&l+XNHZT#8jfV*{
z?twF*(B&dqT{g-JO5t?MFvOc@!o@D)CtSoiO~iP(FvJZrq13hN^*HL?XwZMSm_`Jc
z%;<2{^Ks62Dl*)6&Vtb{?pKlfHT0P`T-<${bCb{D0m^V^SHJ??bN5t0*3i&io)H{l
zJG^}CO7VaVWNqL=Yylh0vf(03Tfz(m9rYZ$I<S{{G*$t0<Z6_hss~g;q8?rWvDCV>
zf|fm%q91(5+QA++5I$#v;R{v{``9$tZ-z{<ow5JR2=FLtMsU0mcH!K%RzbnVp;-4&
zaBf*_8I}}{Wg2FP2G$=}JTt8gP6ChCMg*>7w$86^w%<>?EE`(POM)EJt6`v@=4px@
z!Cl1Nf|@z>2-paz+?7&tSlShlF6K~V`HLrWii;U%z?f*^$=x%!JbbrE><}t8MWQG5
zr9tYov)#pnNbw{t_8eH2YhY$v8J^5|+&31Ftbw7pXjj!hL9fOnz*x_q>97SQrYZu~
z1FP6=RP&)3eGIdEHkY=I>jr%SRgHeV2EF&NK3xN8p*o$l3CaR=<h~@!NKi~P;7wg*
zB=Yz3oJeZ;MEZD6q`xQg1!>{jgFU$i1P%Vm9JkL2clHW-!ev3?b4cv+Id}F75fAh%
z^)kbU<(<t2c`|o)G2@9Eea<_*4ff>j<l;8WGrTOqpZ2t8n{eigN>A=*f+i}su9-b+
zyog~)a1T#>Ub^k+PNX4wUhX-m3&O|c&ZQydp#e={0ie*8h$k46Eb<)ZFc?mIC5*eP
zi`y`}Gjfw^x+^mqZQ|T!7Fp)0X#|WU;`)2%LKMfXv~6Cnh7Ht>!6Pqbg;OvJMhCZ2
z3Ukh}?lxY$^AkKeCMGH?7XO`Q3Gqm<({nl7g_~`==WwndvjHACW?$cI@ZKk(D}&4V
zI1;<I$4*btS9z{Te7H@$k5N~9`ohJoNd>mY3tp=EE>Dx>xVX)7dL6<zjPX?6Bb+(o
zbx+lA1)X94{ATui+lv^+2JM-C1LXB=s^G+X_>;-|EKdBJC5Vq)jKOnyW6hF)D|by$
z!V;8#-n9J&NNPtxK5)n3!7wh!aHzS;{iyP5RQV06Jn+A(ycWhAcOSN(nG+mxaW(A?
zm=GWi+gN?Vo#2Qk^F&YPelF&K6CC&CF7@Q@>*6*}@T({DBv0o4;mjGQJh>;s6k{U)
zjm=Ed;7ME-APzIppfFXqO+)J%vTwO3^I#XVF_Dy>+*3Wdhq$=OL{fP%^LIQqS#j7t
zkr~lm+%OHM8xuXb5N4v|rZ0l7{b_=MVs&g&but!a1Wy)Yf$9`Gb&o}s+F-dH!DLE8
zGM!atL4}cc$09634(mNA(bOn(`eCf7Nn(+aLx<iXBPUSrpq!Y<XsT5zJ@@~m;q&V1
zIjt(Fra_=^szd#xg>E9h4>xHyPwv@GCdKmznVx|R*ZUV2Cj=+O>Bk(G-bS_h{5BEc
z57nfymy0>%MEDbr2Xk(?=DwcH_*;+%bB7RfTT{2o%RIU9_aoxgJ8p7^OF6`ocy3^2
z^wH+T!##)r@Fy=Lxzl2(LU-7-7;=k4<?r|=cNiZ$7Duux%uP+e9q=&Rz>4t{mWk(B
z{6&mJZCC<BP}1c~Agig6w$X{G<R!Q4(Im46%x^je;m0x!%;ufW;2j>vi-5qD_ajSS
zia=qp%iyH((nBnX-OD=R+3_5H#Iul%gd4FX*v)JREQCdPgd4(cMt|Lis!CX$xt-;F
z1G^@`!RNBY0S<m6TM*#j$!vdsgU7IA0S=hmr~~VepbkG=)u;q-gMvyRu~7-0Dr>=y
ziUoYhs^BJUN%)q{fW_Fd;2<l*n~IGr0*<hWa5J`Sc%LnTThJr`JJ^<hVk~7x2kc|g
z*LoJ;z)JXeHk2tpUb%6sY|NQ=sSQgx`fvC1XfIrl<Jd;iHLy#Ncs^5Jw12Q^;nYb1
z(O%SllANjdA8aG74+H#VI^fM>J@)B3Wt`q}8+7x5%-_XF$_IIa{ucZEHJ!Q~y72Y-
z^5sy#C+QoPLt!E>>@12}@^QQPUpKOjd7Wj^PG;+umP2ONIrC<`d6Szbcd`yov1^}W
zSB<G>SGVL7^t;PpP|u0xbrHSph*4t%3fRR<uX0mP!V4yPY_6Ghu?bsQsgXQck6Z!i
ze2T7I0m);^c5$_lWw{y4iBW_Y^WQxC+g6r16lq%WX($QjH(l;#OYlK_H1>C30Lp-3
zDB)8v)>3uy6n=BGe+7*6XUp$+FL{6uk>B$X^1r-F9^?z;A-+WZz_-Z5e7pRSKP8Xw
z7vxWTpFGMB$z%Lyc^u;93CNT`Ll5~o6w5!LM4p5`@)Qh_r(w7}1Eb|xxJI5c`<d>%
z2k{edEBcXO7eOT4hQ0-?5L)~dmct6%OnuA+=1SUpFwW&%c`Ukt722@sEE{P2YIT=_
zF8II4MQ2#<dAdhHzd!(MC~x5B113=Xz7ldK%@~uHmR`?ij=lx>*7FB$;RrvYUl;-w
zFjp-W{%&mni?IZ2gSZ3sFR&uapeffxvaCQ}!ExoA;6UNHpmBV9;VMY7Z-;e8>84eX
z9;&3R!2JYC!As|+a1zxdB`PSH!+NO0lX_#KX$<%y!t#vek$eL`$-x~s^*c>K^@mFt
z84BJec3YkoAN)Db$e$O^?0pZi0d!!au>!qwHN=PR_AKwakh{<%0C(SdbT#0gs6V<|
zZeJ<zV@G{oiq)MycdY+B#Tp*j^tGcdY-=ryY5rNGey*)`wH&wxn0{p~wA5K^AVL?m
zv+~6Db)c))LAow%Z#jBtJ1bL8dj!cwX5$x-y&plcy_r^={`3Wii@}kYe-AEh{UVZ7
zA3_q5V-<=GAAlb70ScMl&T64o?S<H=5bbArVUE>bR-Xd4;Roc{iyX4n5rJM+j8HnS
zJwAFWNub&kk*Gh(v|8!q?X6T<*9M=>JjOHhYe*t4$8X<a^F=a#fo0{%gWd7vhs7d6
zztmks=(}=o0+Vxb0&scYMLG$N*H7a!7+tMD&b4~WxRv;B3)vvA11fEJQLs8*n5^eg
z@12iLFGIhHJu#+1==28EmXv`OXw%V36%{#B7v{(${a9;!RF04Mc)hrzHBg)#EA`7`
zWww4i1BY>{qm?Pfzan+(S7fH{o{hwp<stDgJt@x`An)HRS>5cFiTaleoaxw;GA=q~
vuZK@cJ^iF?OU8`TDPYC9j9I4pbh7%1cYT&#>$5WS!=0?Sh#yku6UF}lsOc%>

delta 10862
zcmb_id0-Vq((kUCm+8DrLI}BUju3KS0tf+&+&2OdbVWolNC;;V0*eR214I-M6i`}O
z#anh=FI-K$*Y#X=6?Io(y#NJ4#rxtCwyNjlW%72u@2^Eu)jd<cs_yFQ>SNxP-O11F
zN{+7@i&jtwCBzfc`)!QxhRxIat?2uBUs=9KU0B_3%%~fW?`m$?Bd<$=@_zc{Px2$?
z-)CZa^G`p^&pIa;EUg&R<-IL`|6-vK{Ff<&&>bI<$@)|VwCh@K!_7fC2X3)pYfvEk
zAy$;atxRrXvW>~@N|$Cr9#UDKl?kgVZFrHNzr^HaCa*Ag)rP+ZMGv@zE3YYdo%tI~
z{-NMa=5G<vP8MW^Y<Mpy2D;YX=U2O!>}K+BCLh@FVNi^41wP^@A2a!c$sSvp%3d4x
z1;u#R-lzQLGaL2?#Y7i%K*8rhI0%Pq_#!CIgIkz<87tm_uT(l5=Gxb>^!^(K-v;43
z_+G&eO23c=)5j1DlNF`}WfG<;OjDSyFhgOc!YqZ^3Ud_ZD$L{e`9bW6ofLLf*hOK1
z!mbJngQ66>DJ)XDE*m<>l2@>g!V0CgW<y<wVHYaANZ~Yv7b~o0@#zXLv2jLFyn{7N
zW-6>@K1<<j=5v_LWyyI8=LfM4>zOoAKX8GKjR6{0{uhkJg}BJZ-v`B~+~-elF}VgW
zRd|`g%N1Us^ffuqKE#?<C|t?>Dut_<U(GgFE4)VGwF=iLTx;Vx8?WPc>x1GLx9lix
zU~;{UHz?c~6sK^L!p#bAWTBfBZc(TeQcNY`Eef|P`~$zZmC0>Nu2PbWZ3~J}Dc;V{
z?y&Jrg?F)djl#PX-lOndh4(4kuJDii_)iM&SNLaz57_u%Q0C!7Ha;Ab9hpB8l%1LH
z2+9I{Na3S#_!oRk;o}OMY-|q75`2Pt|0I*Am^{to8HJ%|6+Wl%uiS>`nf#5(3rt>Q
z@{+=rnZLs1RVEbMrT7|q=yfi=!Q>xI-emF?lbuZ7X7UaNCcew$p9=rQ-Jp0c#rGBN
zvT?VK|DF?+<H()zEc}4Uhb;1u!jBbx!o@ub_bS|{@Kc4KDcrB{fWpre9#nWpA;t4}
z{8Hgp3J)v%TH!YezqRo@h2IBd4Q{gWhoGE=Kic?{jX&##ow8^1vw3iiZM{hD*6p`J
zu^w0m&2mQbuZ55R_6<^G)z&Q1`~HT>&H2SphVqx@s_rl!<(2xma+qrGw&h8((A=>H
z<jL{_@^dNL*Hzcm*Uv4RSTnPxv8L{lnuU%$fWJEOCb`9tTHfr)TbOK>H|vI;P!alr
zyjcosUiHF7j=Yu2!jZQzk6C#6!bLUn9l4F=Z)b7`lRKH*#pG@#_b|CvHahY?x!sX}
zWb!8__cQr3lLweQ$RZDM>0u_1FxerumpbxMF8+ndV@w`r(!`{h$rDVTWb%~M9pkO^
zT1P(3x}IV3ER*NBxryP9kEi}m7<*PdiTzbRF9kNtnCZyBaVuYtjZy?=QdfIB@<o2~
z62o6+@(NpcmC4_^`kLIX-{}RZX|KyS9Qh9>Z!&p{$xivk=I58EWH|C|`HmysW%5rZ
z|6=ly{FqaHYS59N$UTnyR5ot82-7ni`I+1=MR5G!@e?MDoi=RrL>rGdcodIq-jSV^
zfyeQLgTHZCaQf^z2TuYxc*+v`(bJ(SMmiQiTC?g6mYqa_Acb9X*`mhkOBOknWd(F@
zB`k}Jv25E?jumSK_4Z2Wi>u`4`qN5SfN_NV-I?%_M8~pq?}2a+rqH{O2Eys;0Hirq
zJe2G7K~RhdR2V)8%HrFReJha&tS8A8gCI4ry<;U=9ULoJQp1m1DL3teVr5&Yj+JJm
z>mK`{G*%w4$S)tr4|UaLn5nxBhM|Rh0HlaLx2}5L?Ap4T84h&j2Ql;O7tNlz(1Gr3
zD3i=?8VsZ3^QcBjvK%X0e!lt6`RQARKq=->6Uv4_(cs<;HDg#E<_n~dWsa3=*;H6Y
zf2TQ6z|HPrtW0aDZd_R7Kv#ZL2%tB<CsXu8Ltq%@$s_vbA<z>KSoykmC`^v|#j!d<
zn!a`@OiVjqb#mk``Mx9HW3rpczvcV-`=Kx}^n3CdoCyOQh=rg7r?K(|DCb_)G}dw-
z%b?tWcz(H<WaX<=wK`i}99Y0FN*ExSoX*uoCZ$|m$Z`ebU#qKCNS!HltZr731D&`t
zslkag)iY+-)#`_b!GQb&R<UETCaY8)A=&;E6vJDBqC3Sv{BW3}uN)2?@iYpNJBGvc
z)YY=hgI`P0x%D!pdE5xtWQFDkDd&m{M5df4=To_k7A%nUB10~ajgk67sxR`?e^2$r
zQT0nb^~>n}<(~Q#RA1t$UrF_)Tn~i}E~5t_kHB)Oukh4YQvE7Voeiw=)Y-t*p89I4
zU*oB>{I#As%deq&sHMki>A^aW04uo8Q)dP1J#|*F!Bc1X>pgXrzrj;y`Hi0XCaQ0i
zjiHE(ZzQ)<nl|pIOl>$NYr`pD8&2uka7xvNQ>r$clC<HJrwyl6Z8&9X!zp1KPC46f
zO529tD?^Q@K-t@HO5lc54mX_AxZ#w^4X0#oIOTK0DWw}uS>13-?1ocrH=MG&;gsMF
zryOrMrFp}fLg4~sd-H%2zTuSf4Sy<Zh?2gkQ{Fe6QorGp{SBvVZ#boT!ztSvPKn-d
z%Jqg*x;LCMzTqLt`ldjM-*C$P#8*m6{f1NaH=Gs#!)XmLoR$H@X%R4-)&aw5DKMN?
z1H&nO7*4swaLOcxQ!-)RpR=WTKr4sgv}{mJ?4>fuee_2wh<Q%y2k{zuM@v1?gUy@9
zLmbFM&D*MI0}%Rx-lhr=LLQS}%CCt3Nt`VYQ~#y>n%*S}I3?c{V9Rf8`K=J`==pc@
zd)7p6%LJtm;(_AFMBouBkYb`S=HW`o${kn0RZg%U<d0$S05UtAEQ^@plE+10M+6?^
z_K4LYAg(bgO#4JkSOv+;pD5zVviw>8!piiDiI5gz%Pw>+LD&7CP|6<!`c)nYx5L|s
zeh)-i;dWw!Zzqn*V=mYkPbMuo)@5LG+9b$P<k;R9!0p(B`(7+8`wr1Z9+xKsvQSGG
zsPvmWN$miYCQs2&f&5wrroqIJNIJ$o<~;AyJSWP>w6XH>F`zFvAsT$|e}jW$+6p^$
z?=(p39dRnfMS%8ak$CZr2Lz}K<Wz`p4ePvzWJ;(t6ll}t*G${DnK9ALe9{-Je8(@X
zb^JP;0i<2rf18fw6`0+t2Tg|}?NmdSJfyo-L*BV8d^&|s9vz|Mcov9)NU0=J5~3IF
zICzHj5a7z0+2V*mb|4Tq!Sm%?5zA}2F>`A*3?K_1lA$qdfda$Vf)L%qdP~~qjT08>
zji^m%ZT&D(d-e4Hal<9FbCZQWIvt)T3r}4FBTcFpKyzD&hr@B~h`S+?Vq1!G5rXzq
z0#PXzLXwI1YC6()ppqr(#9&CKQouMVPXtbhD?|*Z6cFD3g@}ubkEb}IeWn}zawfq!
z8|;5tFP{O2s4cT=U}Ob5(M=2cBW>ADVMuLh8Hf`iophI#mGviusgPy{;$RIlkdDbS
zq4U_5fe>O^%e)lRAl(dlXHjW-24aNBaQhS13%`kpZkyvsZ${WpugruIq&L163M*M}
zKF#+I!h%f5BJX60_K*!ZRH6ZPZ5eI0fGIZQwkJijg*6RqsD-k^NDhUGp6JaDw^BlW
zt8gcMs21Wqh$+5^j$y=vo#>sIFx?{^i0G${qwD@Q3EqyHm@v~9OdItucomHtcX>3e
z4keTb^Pd*yo!>0<BeP%-IqUmb(4h-E%UvQydmMv87LZeDH{RluWL-HMvON7<>@nWz
z8QR3RK_uuqX18(A6}~Xq=o{F-X2Vo!S&uo;p`f*8V?8aSXe_j-X>=+uh}m<X+wXY1
zl6V<_^+Lf0kq9@46xeJ`aU{f3;^=1LN{F<ums4nhn`b5xVg)49OyR5{p(w2T;2h{m
ztx1>*X(ns1Du+}RhgGG7Rdpa$=2?rXcC_5Ms%WRr&E!D{4@REeHxH64Xh@tUaekAS
z8a=ApJs!#nTabLvF#UaMF7)w?!owcnf`~Ak&Wv!{JgDx=&Pm|q51Pq2Ptfh&Qyx&2
z3o<RIQ7@(L)4|A~ADIV}Jm?n*{gNnvS3Kw;ZP6SvbcQnMGv>n*a$fU%$Q#!>)8~ZS
z*oCendvLZ8qeLg@NhJ_tL_YMQeH@y`P?JxJ$?i;-@w~XvsjEeErarR{s!9Frb^kN<
zbV@bq<LaS_G`~;}Wqr8eZkoF{-0(cHnzj=yyE9rBPKe<&)VW!msD*TBHY#Y4^7Y6D
zNbS{{NsfBl);dxZE`rQdhvTDy<0D1i+W>_#B50HzbVax+68gFxA+&c7mEO7rn<W@y
zJ>b?P(=Y7K@&%A;_HM3ldtc%HzLq*hTT1Z-pWzG6_kgEyp27@Ybfqsk&x1DL9A7XU
z_sxiOX@qe!sEv(KHjoFEr?L>0VNX&C2AVAti2k&>-Ol?#dLI<D=I2Ykp`p<o<nDO%
z-bUzU9F|J^6oGvNV1;n7uSmsyA|LykhPYAFDQo~Jv<I^RqgDw#K{s}>v7~f}Z|^dr
zdp9Jicky&mAoNfdt#4o0W(>zfNxQ<sd`)FXoANGst<f^v7o6h()0vpO;a&b((Ib2}
zi@6@O0sF6pBYnZ$qro};%i*Z-sJ*?Yjf>{^xk&Et?PaHEQ+2-JM}5JaJz&>Gjeckt
z9mxdVxttiaI-qyiZT0UMUl*O`G40L~@9Nu%9ve<uaxtW~t}VFCx1DFfI5VRLkgz8~
z`5!jv291R&R=u67jBeLzUrXa*0)v~MSWG*7n!GPv3L}TLX5{dBnbPJFsDiV@b1OqQ
zCq-WyCQ?$C@d*iWdF|-mQIVWP0)KIvqbFR}CKev_H9yIBKzc_Hi2u>*9A9vI57>=`
z$9+wo>x)kEpv`dlAIyH|8w-7+eV+5IpK1A6*hR}PZ{s4rQ=EkJd`<O_HubWvrB{4o
z;|vd&mQ>0pulb>2ayVuyu7Kjo$o>!C6gIvk5^<+U#&<oy@X5bro0Ckr^Zc-e+uCY)
zpET?u4ZBIhhyTBZDPaw`gtmW&^l3}z5`^Q^-B0_Nya1;9OozeI<9f(rwsn?X=nEdw
z7CcqpVPEt`zUX0X(E`8o1y6IqI{(TxLGhEX@WsBChDTdE;tQ^Z={_DB<pI0%=7cZ$
z5?}Oa51J#2*7DX-YW437U+|b{aE|4Ru5r=&?4{luYAHWq`C53AZD=bT?<Wj1p_W~u
z-&hKBDatCAL9gN5e!jyHVoGESGX-Xa&m~jLW`<AgQ(!jjI@}F)jIc72)A_7B2j&`$
z4=;nEr13}6cm`{9H$^=o8heH}MMCs6dfZJB=}Cx<=e^6k@X}fxYSZUlzTKD)b;js6
z5^mjYS{1&+^^x2_2NE7L0hT8*m{B~#gM@IKC4bQI381cR!Bd4b(1UJ`0s7A4gDz}~
z7S>Q-Fntg*J-mK-n;wqz6{c@SK9-82EsgaBFY?X1Wgf6eoz?^&Gyr`EGwmF)0_th_
zn^r*KfY#x^C^GyPh0n~f;#_mH;|@Qaluw8ObS}%KQ!jlcGpe#y!dOyOyAtvu+iRPm
zhWET2hy5c;`@^M?(-Zwn=7g|i$uZhTOSmk2KmPJcn1Q6Mcs2ZL(&-_QCT<fYa5<e4
zl~^jyhAXIq;wmu~mcW%X4aSPAC_rx}Q-j53bDt~eg8Bj%As34kE<#=@mbeI+F7~(x
znIOJ(5wM`e25lz7HbQV-ixJ3cF#;(qMqsyC1WU;<!pEW>mQl%qPsJPvQ7HfiMGY*c
zQVbjx7sCoF+3=QF1}o|68DNLF(=|+M+$7gOW_j8ql1_-h@{|}SREX~AB<*C(G59JQ
zu97qt`ko?x@dp7alj}~1ek6RkP=2z<#jI#)y6|Sze$xM?IakT!Vk?i20QxTNMl4RE
zI$yRF>u0WoilBVAv`L<A-6LPuzfoOo(!JNv;r?P>y9UZ+xz=l-T#1w2v1bQ4Rm-*X
zvG&j!NYce?p>O3yP4b+rqGz%7xG17Kk*e-i&n7wTyhp`_?q${F4)S8UC*jBDwxmf+
z+b*h&u4=t)Eo6tL(+fhL*(5%u+Z6M@#Js<xNvf?P&#P(%tHM@QrLUeW+%8hakxU0U
zld#m2TI&h1k{%?*Q+)@enjHGX+fLR{8s&Y)eWHk#kT*n*^`<Dd-VwUrIv5tpwLX()
zTKnZ_>wv7XK9}{@LAk^_B$r!X$aU72@_y?p*<>A-&s$&1z1BDKu=On{>qq(wd<5EA
zKS84PGh|r5K#uh*bhM5^k#*e6UA_t*E&mEvQ^=s`EqB3c3L_BZat~Ys*V42}<h#w0
zu!iRuZA08|NPKuL63z)hze-#Cn;2z6XkCvo`mcZQW1{dBe<Gp?AOhF%l0ZKfFcZ1z
zIw%OuI=?t8yIIbjv<l?*Q~Os*qPzI>kn4d=zJs<-3fnp*5-2h{5${Iz(*kJP%*#lu
z$a9wrYF4R91DwRmJk3kws}rD!97Er7a@`cL0j{T3gxu>gG)MF}8rpGLJ+rMDJ)>6Q
zKup9M-CIWzOkBiTT(|`iZh(zsD*6~w6|FF66gs(P!3i%4?k7FiL_XRaUhX!o_bzvV
z1Q+E6Q74`ZF9vUsaO*KCkmM%}H^NN>iEiIH(d|q11#cn8(oZC=8MOiE51U7AupV58
z^dm=oBqPw5f9BZSD>E=5Hu5V+y|!~;X+qndGwPHsfywaAdR^ZoFx;B75rpn>6J0~?
z-3T%Ifr3DLed+zsE{;r@fA_3N&?mbDMp~~vOOg$T=tqzj76i&=rQT2w7$ukLT?K*R
zn0kbMx#3`nKD}$895Wr{T4z@0=q+6Xxi~ou_0BZR)=qXHQGeDokSXb>mt6x@@I#5N
zDhv#@*7rrBcTB=${pZ3!x%J1fZZ%sUrz#$;N8Qjs2AXpSKdD<F6PniQ`fk*PyPiYJ
z(&g#;?QYbr1zHx7)>Fszq@uu3*+FY!mNnQCo8P`8MgLety6>JszuUR^N134q6%)kL
z3yK3{<W2gS;=l-N#0Mn0U|ovNDxr60>9a}#6J(sewInc1HtKyPfwGwEGx(Q-{{g89
B__P23

diff --git a/src/nanook/NanoOK.java b/src/nanook/NanoOK.java
index 6f36e77..d35b695 100644
--- a/src/nanook/NanoOK.java
+++ b/src/nanook/NanoOK.java
@@ -25,7 +25,7 @@
  * @author Richard Leggett
  */
 public class NanoOK {
-    public final static String VERSION_STRING = "v1.09";
+    public final static String VERSION_STRING = "v1.10";
     public final static long SERIAL_VERSION = 3L;
     
     /**
diff --git a/src/nanook/NanoOKOptions.java b/src/nanook/NanoOKOptions.java
index a2259f9..45fce12 100644
--- a/src/nanook/NanoOKOptions.java
+++ b/src/nanook/NanoOKOptions.java
@@ -1246,6 +1246,9 @@ void readProcessFile() {
                         } else if (tokens[0].compareToIgnoreCase("Reference") == 0) {
                             referenceFile = tokens[1];
                             System.out.println("  Reference "+tokens[1]);
+                        } else if (tokens[0].compareToIgnoreCase("Sample") == 0) {
+                            sampleDirectory = tokens[1];
+                            System.out.println("  Sample "+tokens[1]);
                         } else if (tokens[0].compareToIgnoreCase("Analysis") == 0) {
                             parsingReads = true;
                             System.out.println("  Analysis "+tokens[1]);