Browse code

Add low-level ..gsnap() interface to gsnap that accepts an argument string (instead of building that string from its formal arguments, like .gsnap). Also checks for whether gsnap succeeds.

git-svn-id: file:///home/git/hedgehog.fhcrc.org/bioconductor/trunk/madman/Rpacks/gmapR@73839 bc3139a8-67e5-0310-9ffc-ced21a209358

Michael Lawrence authored on 28/02/2013 23:03:33
Showing 3 changed files

... ...
@@ -10,7 +10,7 @@ Description: GSNAP and GMAP are a pair of tools to align short-read
10 10
     to work with GMAP and GSNAP from within R. In addition, it provides 
11 11
     methods to tally alignment results on a per-nucleotide basis using 
12 12
     the bam_tally tool.
13
-Version: 1.1.9
13
+Version: 1.1.10
14 14
 Depends: R (>= 2.15.0), methods, GenomicRanges
15 15
 Imports: IRanges, Rsamtools (>= 1.7.4), rtracklayer (>= 1.17.15), GenomicRanges,
16 16
          GenomicFeatures, Biostrings, VariantAnnotation, tools, Biobase
... ...
@@ -51,7 +51,8 @@ setMethod("gsnap", c("character", "characterORNULL", "GsnapParam"),
51 51
 ### Low-level interface
52 52
 ###
53 53
 
54
-.gsnap <- function(db, dir = NULL, part = NULL, input_buffer_size = 1000L,
54
+.gsnap <- function(db = NULL, dir = NULL, part = NULL,
55
+                   input_buffer_size = 1000L,
55 56
                    barcode_length = 0L,
56 57
                    orientation = c("FR", "RF", "FF"),
57 58
                    fastq_id_start = NULL, fastq_id_end = NULL,
... ...
@@ -107,7 +108,6 @@ setMethod("gsnap", c("character", "characterORNULL", "GsnapParam"),
107 108
   if (!is.null(problems))
108 109
     stop("validation failed:\n  ", paste(problems, collapse = "\n  "))  
109 110
 
110
-  ##TODO: figure out how to do this programmatically
111 111
   orientation <- match.arg(orientation)
112 112
   batch <- match.arg(batch)
113 113
   mode <- match.arg(mode)
... ...
@@ -116,7 +116,26 @@ setMethod("gsnap", c("character", "characterORNULL", "GsnapParam"),
116 116
   
117 117
 ### TODO: if input_a is NULL, or split_output and .redirect are NULL:
118 118
 ###       return a pipe()
119
-  .system(commandLine("gsnap"), intern = version)
119
+  .system_gsnap(commandLine("gsnap"))
120
+}
121
+
122
+..gsnap <- function(args, path = NULL) {
123
+  .system_gsnap(.commandLine("gsnap", args, path))
124
+}
125
+
126
+.system_gsnap <- function(command) {
127
+  command <- paste(command, "2>&1")
128
+  output <- .system(command, intern = TRUE)
129
+  if (!gsnapSucceeded(command, output))
130
+    stop("Execution of gsnap failed, command-line: '", command,
131
+         "'; last output line: '", tail(output, 1), "'")
132
+  else output
133
+}
134
+
135
+gsnapSucceeded <- function(command, output) {
136
+  if (!grepl("--version", command))
137
+    any(grepl("^Processed \\d+ queries", output))
138
+  else TRUE
120 139
 }
121 140
 
122 141
 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
... ...
@@ -198,10 +217,11 @@ setValidity("GsnapParam", .valid_GsnapParam)
198 217
 ###
199 218
 
200 219
 gsnapVersion <- function() {
201
-  version_text <- sub("GSNAP version (.*?) ", "\\1", .gsnap(version = TRUE))
220
+  output <- .gsnap(version = TRUE)
221
+  version_text <- sub("GSNAP version (.*?) .*", "\\1", output[1])
202 222
   parseGsnapVersion(version_text)
203 223
 }
204 224
 
205 225
 parseGsnapVersion <- function(x) {
206
-  as.POSIXlt(x, format = "%m-%d-%Y")
226
+  as.POSIXlt(x, format = "%Y-%m-%d")
207 227
 }
... ...
@@ -31,7 +31,7 @@ getDefaultGmapPath <- function() {
31 31
 
32 32
 ## gathers command line from non-NULL args to parent frame
33 33
 commandLine <- function(binary = "gsnap",
34
-                        path = getOption("gmap.path", getDefaultGmapPath()))
34
+                        path = NULL)
35 35
 {
36 36
   ##get values of arguments of function that called this function
37 37
   parentFormals <- formals(sys.function(sys.parent()))
... ...
@@ -94,13 +94,18 @@ commandLine <- function(binary = "gsnap",
94 94
   ##long-form args separate arg name and value with "=". Not strictly
95 95
   ##necessary since getopts handles them, but has caused confusion so
96 96
   ##adding
97
-  namedUserArgs[nchar(dashes) > 1] <- sub(" ", "=", namedUserArgs[nchar(dashes) > 1])
98
-  if (!is.null(path))
99
-    binary <- file.path(path, binary)
100
-  paste(binary, paste(c(namedUserArgs, unnamedUserArgs), collapse = " "))
97
+  namedUserArgs[nchar(dashes) > 1] <-
98
+    sub(" ", "=", namedUserArgs[nchar(dashes) > 1])
99
+  .commandLine(binary, namedUserArgs, unnamedUserArgs, path = path)
100
+}
101
+
102
+.commandLine <- function(binary, ..., path = NULL) {
103
+  if (is.null(path))
104
+    path <- getOption("gmap.path", getDefaultGmapPath())
105
+  arg.string <- paste(c(...), collapse = " ")
106
+  paste(file.path(path, binary), arg.string)
101 107
 }
102 108
 
103
-## at some point, mxbay want to customize this
104 109
 .system <- function(...) {
105 110
 
106 111
   if (is.null(getOption("systemCallMode"))) {