diff scripts/dropletutils.Rscript @ 2:a8aa294401be draft

"planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/dropletutils/ commit 4d89eb1eb951ef094d1f77c46824d9c38be4445b"

--- a/scripts/dropletutils.Rscript	Mon Aug 26 05:06:39 2019 -0400
+++ b/scripts/dropletutils.Rscript	Fri Sep 06 10:56:16 2019 -0400
@@ -58,35 +58,44 @@
 
     eparams$... <- NULL ## hack
     eparams$m = Matrix(counts(sce), sparse=TRUE)
-    
+
     e.out <- do.call(emptyDrops, c(eparams))
-    
+
     bar.names <- colnames(sce)
-    if (length(bar.names) != nrow(e.out)){    
+    if (length(bar.names) != nrow(e.out)){
         stop("Length of barcodes and output metrics don't match.")
     }
     e.out <- cbind(bar.names, e.out)
     e.out$is.Cell <- e.out$FDR <= fdr_threshold
     e.out$is.CellAndLimited <- e.out$is.Cell & e.out$Limited
-    
-    # Write to table
-    writeTSV(files$table, e.out)
-    
-    # Print to log
-    print(table(Limited=e.out$Limited, Significant=e.out$is.Cell))
-    
-    # Write to Plot
+
+    ## Write to Plot
+    e.out$is.Cell[is.na(e.out$is.Cell)] <- FALSE
+    xlim.dat <- e.out[complete.cases(e.out),]$Total
+
+    ## Write to table
+    writeTSV(files$table, e.out[complete.cases(e.out),])
+
     png(files$plot)
     plot(e.out$Total, -e.out$LogProb, col=ifelse(e.out$is.Cell, "red", "black"),
-    xlab="Total UMI count", ylab="-Log Probability")
+         xlab="Total UMI count", ylab="-Log Probability",
+         xlim=c(min(xlim.dat),max(xlim.dat)))
     dev.off()
-    
-    # Filtered
-    called <- e.out$is.CellAndLimited
+
+    ## Filtered
+    called <- NULL
+    if (fdr_threshold != 0){
+        called <- e.out$is.CellAndLimited
+    } else {
+        called <- e.out$is.Cell
+    }
     called[is.na(called)] <- FALSE    # replace NA's with FALSE
     sce.filtered <- sce[,called]
-    
+
     writeOut(counts(sce.filtered), files$out, out.type)
+
+    message(paste("Cells:", sum(na.omit(e.out$is.Cell))))
+    message(paste("Cells and Limited:", sum(na.omit(e.out$is.CellAndLimited))))
 }
 
 
@@ -95,12 +104,13 @@
 
     dparams$m = counts(sce)
     called <- do.call(defaultDrops, c(dparams))
-    print(table(called))
-        
+
     # Filtered
     sce.filtered <- sce[,called]
-    
+
     writeOut(Matrix(counts(sce.filtered),sparse=TRUE), files$out, out.type)
+
+    message(paste("Cells:", sum(called)))
 }
 
 
@@ -111,7 +121,7 @@
     bparams$m = counts(sce)
 
     br.out <- do.call(barcodeRanks, c(bparams))
-    
+
     png(files$plot)
     plot(br.out$rank, br.out$total, log="xy", xlab="(log) Rank", ylab="(log) Total Number of Barcodes")
     o <- order(br.out$rank)
@@ -121,7 +131,7 @@
     abline(h=br.out$inflection, col="forestgreen", lty=2)
     legend("bottomleft", lty=2, col=c("dodgerblue", "forestgreen"), legend=c("knee", "inflection"))
     dev.off()
-    
+
     print(paste("knee =", br.out$knee, ", inflection = ", br.out$inflection))
 }