diff --git a/R/plot_residualError_OtherCovariate.R b/R/plot_residualError_OtherCovariate.R index a85149061..ba187036e 100644 --- a/R/plot_residualError_OtherCovariate.R +++ b/R/plot_residualError_OtherCovariate.R @@ -54,8 +54,12 @@ pdf(outfile, height=7, width=7) hst=subset(data.frame(e$Covariate, e$nBases), e.nBases != 0) hst2=subset(data.frame(f$Covariate, f$nBases), f.nBases != 0) if( is.numeric(c$Covariate) ) { - plot(hst$e.Covariate, hst$e.nBases, type="h", lwd=2, main=paste(covariateName,"histogram"), xlab=covariateName, ylab="Count",yaxt="n",xlim=c(min(c$Covariate),max(c$Covariate))) - points(hst2$f.Covariate, hst2$f.nBases, type="h", lwd=2, col="cornflowerblue") + if( length(hst$e.Covariate) == 0 ) { + plot(hst2$f.Covariate, hst2$f.nBases, type="h", lwd=2, col="cornflowerblue", main=paste(covariateName,"histogram"), xlab=covariateName, ylab="Count",yaxt="n",xlim=c(min(c$Covariate),max(c$Covariate))) + } else { + plot(hst$e.Covariate, hst$e.nBases, type="h", lwd=2, main=paste(covariateName,"histogram"), xlab=covariateName, ylab="Count",yaxt="n",xlim=c(min(c$Covariate),max(c$Covariate))) + points(hst2$f.Covariate, hst2$f.nBases, type="h", lwd=2, col="cornflowerblue") + } axis(2,axTicks(2), format(axTicks(2), scientific=F)) } else { # Dinuc (and other non-numeric covariates) are different to make their plots look nice hst=subset(data.frame(c$Covariate, c$nBases), c.nBases != 0) diff --git a/R/plot_residualError_QualityScoreCovariate.R b/R/plot_residualError_QualityScoreCovariate.R index 8dfaf0b3d..7dc885205 100644 --- a/R/plot_residualError_QualityScoreCovariate.R +++ b/R/plot_residualError_QualityScoreCovariate.R @@ -21,7 +21,7 @@ e <- rbind(d.good, d.1000, d.10000) rmseGood = sqrt(sum((d.good$Qempirical-d.good$Qreported)^2 * d.good$nBases) / sum(d.good$nBases) ) rmseAll = sqrt(sum((e$Qempirical-e$Qreported)^2 * e$nBases) / sum(e$nBases) ) theTitle = paste("RMSE_good = ", round(rmseGood,digits=3), ", RMSE_all = ", round(rmseAll,digits=3)) -if (length(t$nBases) == length(d.good$nBases) ) { +if( length(t$nBases) - length(f$nBases) == length(d.good$nBases) ) { theTitle = paste("RMSE = ", round(rmseAll,digits=3)); } plot(d.good$Qreported, d.good$Qempirical, type="p", col="blue", main=theTitle, xlim=c(0,40), ylim=c(0,40), pch=16, xlab="Reported quality score", ylab="Empirical quality score")