--- a/goseq.r	Thu Feb 25 08:49:20 2016 -0500
+++ b/goseq.r	Fri Feb 26 08:49:59 2016 -0500
@@ -46,35 +46,34 @@
   length_table = read.delim(length_file, header=TRUE, sep="\t", check.names=FALSE, row.names=TRUE)
   gene_lengths = length_table[names(genes),]$length
   } else {
-  gene_lengths = getlength(names(genes),genome, gene_id)
+  gene_lengths = getlength(names(genes), genome, gene_id)
   }
 
 # Estimate PWF
 
-print(gene_lengths)
-
 pdf(length_bias_plot)
 pwf=nullp(genes, genome, gene_id, gene_lengths)
 dev.off()
-# Null dstribution wallenius
+# wallenius approximation of p-values
 GO.wall=goseq(pwf, genome, gene_id)
 
 GO.nobias=goseq(pwf, genome, gene_id, method="Hypergeometric")
 
-# Sampling dsitribution
-GO.samp=goseq(pwf,genome, gene_id, method="Sampling",repcnt=repcnt)
-
-# Compare sampling with wallenius
-pdf(sample_vs_wallenius_plot)
-plot(log10(GO.wall[,2]), log10(GO.samp[match(GO.samp[,1],GO.wall[,1]),2]),
+# Sampling distribution
+if (repcnt > 0) {
+  GO.samp=goseq(pwf,genome, gene_id, method="Sampling", repcnt=repcnt)
+  # Compare sampling with wallenius
+  pdf(sample_vs_wallenius_plot)
+  plot(log10(GO.wall[,2]), log10(GO.samp[match(GO.samp[,1],GO.wall[,1]),2]),
      xlab="log10(Wallenius p-values)",ylab="log10(Sampling p-values)",
      xlim=c(-3,0))
-abline(0,1,col=3,lty=2)
-dev.off()
+     abline(0,1,col=3,lty=2)
+  dev.off()
+  write.table(GO.samp, sampling_tab, sep="\t", row.names = FALSE, quote = FALSE)
+}
 
 
 write.table(GO.wall, wallenius_tab, sep="\t", row.names = FALSE, quote = FALSE)
-write.table(GO.samp, sampling_tab, sep="\t", row.names = FALSE, quote = FALSE)
 write.table(GO.nobias, nobias_tab, sep="\t", row.names = FALSE, quote = FALSE)
 
 sessionInfo()