@@ -50,8 +50,8 @@ setGeneric("genomeName", function(x) standardGeneric("genomeName"))

 setMethod("genomeName", "character", function(x) x)

 setMethod("genomeName", "BSgenome", function(x) providerVersion(x))

-setMethods("genomeName", list("RTLFile", "RsamtoolsFile"),

-           function(x) file_path_sans_ext(basename(path(x)), TRUE))

+setMethod("genomeName", "BiocFile",

+          function(x) file_path_sans_ext(basename(path(x)), TRUE))

 setMethod("genomeName", "ANY", function(x) {

     if (hasMethod("seqinfo", class(x))) {

         ans <- unique(genome(x))

@@ -121,16 +121,7 @@ setReplaceMethod("referenceSequence",

                  signature(x = "GmapGenome", value = "ANY"),

                  function(x, name, ..., value)

                  {

-                   gmap_db_tmp_dir <- file.path(tempdir(), "gmap_db_tmp_dir")

-                   dir.create(gmap_db_tmp_dir, recursive=TRUE)

-                   cur_wd <- getwd()

-                   on.exit({unlink(gmap_db_tmp_dir, recursive=TRUE)

-                            setwd(cur_wd)})

-                   setwd(gmap_db_tmp_dir)

-

-                   db <- gmap_build(value, x, ...)

-                   

-                   db

+                   gmap_build(value, x, ...)

                  })

 setGeneric("snps<-", function(x, name, ..., value) standardGeneric("snps<-"))

@@ -178,7 +178,10 @@ setReplaceMethod("spliceSites", c("GmapGenome", "TxDb"),

 setMethod("getSeq", "GmapGenome", function(x, which = seqinfo(x)) {

   if (!.gmapGenomeCreated(x))

-    stop("Genome index does not exist")

+      stop("Genome index does not exist")

+  if (is.character(which)) {

+      which <- seqinfo(x)[which]

+  }

   which <- as(which, "GRanges")

   merge(seqinfo(x), seqinfo(which)) # for the checks

   ans <- .Call(R_Genome_getSeq, path(directory(x)), genome(x),

@@ -52,6 +52,17 @@ setMethod("genomeName", "character", function(x) x)

 setMethod("genomeName", "BSgenome", function(x) providerVersion(x))

 setMethods("genomeName", list("RTLFile", "RsamtoolsFile"),

            function(x) file_path_sans_ext(basename(path(x)), TRUE))

+setMethod("genomeName", "ANY", function(x) {

+    if (hasMethod("seqinfo", class(x))) {

+        ans <- unique(genome(x))

+        if (length(ans) > 1L) {

+            stop("genome is ambiguous")

+        }

+        ans

+    } else {

+        stop("cannot derive a genome name")

+    }

+})

 file_path_is_absolute <- function(x) {

   ## hack that is unlikely to work on e.g. Windows

@@ -50,8 +50,8 @@ setGeneric("genomeName", function(x) standardGeneric("genomeName"))

 setMethod("genomeName", "character", function(x) x)

 setMethod("genomeName", "BSgenome", function(x) providerVersion(x))

-setMethod("genomeName", "RTLFile",

-          function(x) file_path_sans_ext(basename(path(x)), TRUE))

+setMethods("genomeName", list("RTLFile", "RsamtoolsFile"),

+           function(x) file_path_sans_ext(basename(path(x)), TRUE))

 file_path_is_absolute <- function(x) {

   ## hack that is unlikely to work on e.g. Windows

@@ -169,6 +169,7 @@ setMethod("getSeq", "GmapGenome", function(x, which = seqinfo(x)) {

   if (!.gmapGenomeCreated(x))

     stop("Genome index does not exist")

   which <- as(which, "GRanges")

+  merge(seqinfo(x), seqinfo(which)) # for the checks

   ans <- .Call(R_Genome_getSeq, path(directory(x)), genome(x),

                as.character(seqnames(which)), start(which), width(which),

                as.character(strand(which)))

@@ -50,7 +50,7 @@ setGeneric("genomeName", function(x) standardGeneric("genomeName"))

 setMethod("genomeName", "character", function(x) x)

 setMethod("genomeName", "BSgenome", function(x) providerVersion(x))

-setMethod("genomeName", "FastaFile",

+setMethod("genomeName", "RTLFile",

           function(x) file_path_sans_ext(basename(path(x)), TRUE))

 file_path_is_absolute <- function(x) {

@@ -154,6 +154,7 @@ setReplaceMethod("spliceSites", c("GmapGenome", "GRangesList"),

                    x

                  })

+#setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

 setReplaceMethod("spliceSites", c("GmapGenome", "TxDb"),

                  function(x, name, value) {

                    spliceSites(x, name) <- exonsBy(value)

@@ -154,7 +154,7 @@ setReplaceMethod("spliceSites", c("GmapGenome", "GRangesList"),

                    x

                  })

-setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

+setReplaceMethod("spliceSites", c("GmapGenome", "TxDb"),

                  function(x, name, value) {

                    spliceSites(x, name) <- exonsBy(value)

                    x

@@ -32,10 +32,13 @@ mapsDirectory <- function(x) {

 setMethod("seqinfo", "GmapGenome", function(x) {

   if (!.gmapGenomeCreated(x))

     stop("GmapGenome index '", genome(x), "' does not exist")

-  tab <- read.table(.get_genome(path(directory(x)), genome(x),

-                                chromosomes = TRUE),

-                    colClasses = c("character", "NULL", "integer", "character"),

-                    fill = TRUE)

+  suppressWarnings({ # warning when colClasses is too long, even when fill=TRUE!

+    tab <- read.table(.get_genome(path(directory(x)), genome(x),

+                                  chromosomes = TRUE),

+                      colClasses = c("character", "NULL", "integer",

+                        "character"),

+                      fill = TRUE)

+  })

   Seqinfo(tab[,1], tab[,2], nzchar(tab[,3]), genome = genome(x))

 })

@@ -34,8 +34,9 @@ setMethod("seqinfo", "GmapGenome", function(x) {

     stop("GmapGenome index '", genome(x), "' does not exist")

   tab <- read.table(.get_genome(path(directory(x)), genome(x),

                                 chromosomes = TRUE),

-                    colClasses = c("character", "NULL", "integer"))

-  Seqinfo(tab[,1], tab[,2], genome = genome(x))

+                    colClasses = c("character", "NULL", "integer", "character"),

+                    fill = TRUE)

+  Seqinfo(tab[,1], tab[,2], nzchar(tab[,3]), genome = genome(x))

 })

 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

@@ -30,6 +30,8 @@ mapsDirectory <- function(x) {

 }

 setMethod("seqinfo", "GmapGenome", function(x) {

+  if (!.gmapGenomeCreated(x))

+    stop("GmapGenome index '", genome(x), "' does not exist")

   tab <- read.table(.get_genome(path(directory(x)), genome(x),

                                 chromosomes = TRUE),

                     colClasses = c("character", "NULL", "integer"))

@@ -189,6 +189,8 @@ setMethod("show", "GmapGenome", function(object) {

       path(directory(object)), "\n")

 })

+setMethod("showAsCell", "GmapGenome", function(object) genome(object))

+

 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 ### Utilities

 ###

@@ -47,12 +47,30 @@ setMethod("genomeName", "BSgenome", function(x) providerVersion(x))

 setMethod("genomeName", "FastaFile",

           function(x) file_path_sans_ext(basename(path(x)), TRUE))

+file_path_is_absolute <- function(x) {

+  ## hack that is unlikely to work on e.g. Windows

+  identical(substring(x, 1, 1), .Platform$file.sep)

+}

+

+file_path_is_dir <- function(x) {

+  isTRUE(file.info(x)[,"isdir"])

+}

+

 GmapGenome <- function(genome,

                        directory = GmapGenomeDirectory(create = create), 

                        name = genomeName(genome), create = FALSE, ...)

 {

   if (!isTRUEorFALSE(create))

     stop("'create' must be TRUE or FALSE")

+  if (isSingleString(genome) && file_path_is_dir(genome)) {

+    genome <- path.expand(genome)

+    if (file_path_is_absolute(genome)) {

+      if (!missing(directory))

+        stop("'directory' should be missing when 'genome' is an absolute path")

+      directory <- dirname(genome)

+      genome <- basename(genome)

+    }

+  }

   if (isSingleString(directory))

     directory <- GmapGenomeDirectory(directory, create = create)

   if (is(genome, "DNAStringSet")) {

@@ -141,6 +141,8 @@ setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

 ###

 setMethod("getSeq", "GmapGenome", function(x, which = seqinfo(x)) {

+  if (!.gmapGenomeCreated(x))

+    stop("Genome index does not exist")

   which <- as(which, "GRanges")

   ans <- .Call(R_Genome_getSeq, path(directory(x)), genome(x),

                as.character(seqnames(which)), start(which), width(which),

@@ -179,3 +181,20 @@ setMethod("show", "GmapGenome", function(object) {

   db

 }

+

+.gmapGenomeCreated <- function(genome) {

+  ##existance means the GENOME_NAME.chromosome exists

+  

+  d <- path(directory(genome))

+  if (!file.exists(d))

+    return(FALSE)

+

+  chromosome.file <- paste(genome(genome), "chromosome", sep=".")

+  possibleLoc1 <- file.path(d, chromosome.file)

+  possibleLoc2 <- file.path(d, genome(genome), chromosome.file)

+  if (!(file.exists(possibleLoc1) || file.exists(possibleLoc2)))

+    return(FALSE)

+

+  return(TRUE)

+}

+

@@ -145,8 +145,8 @@ setMethod("getSeq", "GmapGenome", function(x, which = seqinfo(x)) {

   ans <- .Call(R_Genome_getSeq, path(directory(x)), genome(x),

                as.character(seqnames(which)), start(which), width(which),

                as.character(strand(which)))

-  if (!is.null(names(x)))

-    names(ans) <- names(x)

+  if (!is.null(names(which)))

+    names(ans) <- names(which)

   ans

 })

@@ -142,9 +142,12 @@ setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

 setMethod("getSeq", "GmapGenome", function(x, which = seqinfo(x)) {

   which <- as(which, "GRanges")

-  .Call(R_Genome_getSeq, path(directory(x)), genome(x),

-        as.character(seqnames(which)), start(which), width(which),

-        as.character(strand(which)))

+  ans <- .Call(R_Genome_getSeq, path(directory(x)), genome(x),

+               as.character(seqnames(which)), start(which), width(which),

+               as.character(strand(which)))

+  if (!is.null(names(x)))

+    names(ans) <- names(x)

+  ans

 })

 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

@@ -153,6 +156,10 @@ setMethod("getSeq", "GmapGenome", function(x, which = seqinfo(x)) {

 setAs("ANY", "GmapGenome", function(from) GmapGenome(from))

+setAs("GmapGenome", "DNAStringSet", function(from) {

+  DNAStringSet(getSeq(from))

+})

+

 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 ### Show

 ###

@@ -140,7 +140,7 @@ setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

 ### Sequence access

 ###

-setMethod("getSeq", "GmapGenome", function(x, which) {

+setMethod("getSeq", "GmapGenome", function(x, which = seqinfo(x)) {

   which <- as(which, "GRanges")

   .Call(R_Genome_getSeq, path(directory(x)), genome(x),

         as.character(seqnames(which)), start(which), width(which),

@@ -136,6 +136,17 @@ setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

                    x

                  })

+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+### Sequence access

+###

+

+setMethod("getSeq", "GmapGenome", function(x, which) {

+  which <- as(which, "GRanges")

+  .Call(R_Genome_getSeq, path(directory(x)), genome(x),

+        as.character(seqnames(which)), start(which), width(which),

+        as.character(strand(which)))

+})

+

 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 ### Coerce

 ###

@@ -58,7 +58,12 @@ GmapGenome <- function(genome,

   if (is(genome, "DNAStringSet")) {

     if (missing(name))

       stop("If the genome argument is a DNAStringSet object",

-           "the name argument must be provided") 

+           "the name argument must be provided")

+    if (is.null(names(genome))) {

+      stop("If the genome is provided as a DNAStringSet, ",

+           "the genome needs to have names. ",

+           "E.g., \"names(genome) <- someSeqNames")

+    }

   }

   if (!isSingleString(name))

     stop("'name' must be a single, non-NA string")

@@ -109,7 +109,7 @@ setGeneric("spliceSites<-",

 setReplaceMethod("spliceSites", c("GmapGenome", "GRangesList"),

                  function(x, name, value) {

                    exonsFlat <- unlist(value, use.names=FALSE)

-                   exonsPart <- PartitioningByWidth(exons)

+                   exonsPart <- PartitioningByWidth(value)

                    exonsHead <- exonsFlat[-end(exonsPart)]

                    donors <- flank(exonsHead, 1L, start = FALSE)

                    exonsTail <- exonsFlat[-start(exonsPart)]

@@ -118,7 +118,7 @@ setReplaceMethod("spliceSites", c("GmapGenome", "GRangesList"),

                               resize(acceptors, 2L, fix = "start"))

                    names(sites) <- values(sites)$exon_id

                    info <- rep(c("donor", "acceptor"), each = length(donors))

-                   intronWidths <- abs(acceptors - donors) + 1L

+                   intronWidths <- abs(start(acceptors) - start(donors)) + 1L

                    info <- paste(info, intronWidths)

                    values(sites) <- DataFrame(info)

                    iit_store(sites, file.path(mapsDirectory(x), name))

@@ -108,7 +108,6 @@ setGeneric("spliceSites<-",

 setReplaceMethod("spliceSites", c("GmapGenome", "GRangesList"),

                  function(x, name, value) {

-                   browser()

                    exonsFlat <- unlist(value, use.names=FALSE)

                    exonsPart <- PartitioningByWidth(exons)

                    exonsHead <- exonsFlat[-end(exonsPart)]

@@ -106,25 +106,32 @@ setReplaceMethod("snps", c("GmapGenome", "ANY"),

 setGeneric("spliceSites<-",

            function(x, ..., value) standardGeneric("spliceSites<-"))

-setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

+setReplaceMethod("spliceSites", c("GmapGenome", "GRangesList"),

                  function(x, name, value) {

-                   exons <- exonsBy(value)

-                   exonsFlat <- unlist(exons, use.names=FALSE)

+                   browser()

+                   exonsFlat <- unlist(value, use.names=FALSE)

                    exonsPart <- PartitioningByWidth(exons)

                    exonsHead <- exonsFlat[-end(exonsPart)]

-                   donors <- flank(exonsHead, 1L, start = FALSE, both = TRUE)

+                   donors <- flank(exonsHead, 1L, start = FALSE)

                    exonsTail <- exonsFlat[-start(exonsPart)]

-                   acceptors <- flank(exonsTail, 1L, start = TRUE, both = TRUE)

-                   sites <- c(donors, acceptors)

+                   acceptors <- flank(exonsTail, 1L, start = TRUE)

+                   sites <- c(resize(donors, 2L, fix = "end"),

+                              resize(acceptors, 2L, fix = "start"))

                    names(sites) <- values(sites)$exon_id

                    info <- rep(c("donor", "acceptor"), each = length(donors))

-                   info <- paste(info, unlist(width(intronsByTranscript(value)),

-                                              use.names = FALSE))

+                   intronWidths <- abs(acceptors - donors) + 1L

+                   info <- paste(info, intronWidths)

                    values(sites) <- DataFrame(info)

                    iit_store(sites, file.path(mapsDirectory(x), name))

                    x

                  })

+setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

+                 function(x, name, value) {

+                   spliceSites(x, name) <- exonsBy(value)

+                   x

+                 })

+

 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 ### Coerce

 ###

@@ -93,7 +93,7 @@ setReplaceMethod("referenceSequence",

                    db

                  })

-setGeneric("snps<-", function(x, ..., value) standardGeneric("snps<-"))

+setGeneric("snps<-", function(x, name, ..., value) standardGeneric("snps<-"))

 setReplaceMethod("snps", c("GmapGenome", "ANY"),

                  function(x, name, ..., value) {

@@ -96,7 +96,7 @@ setReplaceMethod("referenceSequence",

 setGeneric("snps<-", function(x, ..., value) standardGeneric("snps<-"))

 setReplaceMethod("snps", c("GmapGenome", "ANY"),

-                 function(x, ..., value) {

+                 function(x, name, ..., value) {

                    snpDir <- GmapSnpDirectory(x)

                    snps(snpDir, name = name, genome = x,

                         iitPath = mapsDirectory(x), ...) <- value

@@ -125,6 +125,12 @@ setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

                    x

                  })

+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+### Coerce

+###

+

+setAs("ANY", "GmapGenome", function(from) GmapGenome(from))

+

 ### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

 ### Show

 ###

new file mode 100644

@@ -0,0 +1,146 @@

+### =========================================================================

+### GmapGenome class

+### -------------------------------------------------------------------------

+###

+### Database of reference sequence used by the GMAP suite.

+###

+

+setClass("GmapGenome", representation(name = "character",

+                                      directory = "GmapGenomeDirectory"))

+

+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+### Accessors

+###

+

+directory <- function(x) {

+  retVal <- NULL

+  if (!is.null(x)) {

+    retVal <- x@directory ## 'dir' is already taken

+  }

+  return(retVal)

+}

+

+setMethod("genome", "GmapGenome", function(x) x@name)

+

+setMethod("path", "GmapGenome",

+          function(object) file.path(path(directory(object)), genome(object)))

+

+mapsDirectory <- function(x) {

+  file.path(path(x), paste(genome(x), "maps", sep = "."))

+}

+

+setMethod("seqinfo", "GmapGenome", function(x) {

+  tab <- read.table(.get_genome(path(directory(x)), genome(x),

+                                chromosomes = TRUE),

+                    colClasses = c("character", "NULL", "integer"))

+  Seqinfo(tab[,1], tab[,2], genome = genome(x))

+})

+

+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+### Constructor

+###

+

+setGeneric("genomeName", function(x) standardGeneric("genomeName"))

+

+setMethod("genomeName", "character", function(x) x)

+setMethod("genomeName", "BSgenome", function(x) providerVersion(x))

+setMethod("genomeName", "FastaFile",

+          function(x) file_path_sans_ext(basename(path(x)), TRUE))

+

+GmapGenome <- function(genome,

+                       directory = GmapGenomeDirectory(create = create), 

+                       name = genomeName(genome), create = FALSE, ...)

+{

+  if (!isTRUEorFALSE(create))

+    stop("'create' must be TRUE or FALSE")

+  if (isSingleString(directory))

+    directory <- GmapGenomeDirectory(directory, create = create)

+  if (is(genome, "DNAStringSet")) {

+    if (missing(name))

+      stop("If the genome argument is a DNAStringSet object",

+           "the name argument must be provided") 

+  }

+  if (!isSingleString(name))

+    stop("'name' must be a single, non-NA string")

+  if (!is(directory, "GmapGenomeDirectory"))

+    stop("'directory' must be a GmapGenomeDirectory object or path to one")

+  db <- new("GmapGenome", name = name, directory = directory)

+  if (create) {

+    if (name %in% genome(directory))

+      message("NOTE: genome '", name, "' already exists, not overwriting")

+    else referenceSequence(db, ...) <- genome

+  }

+  db

+}

+

+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+### Building

+###

+

+setReplaceMethod("referenceSequence",

+                 signature(x = "GmapGenome", value = "ANY"),

+                 function(x, name, ..., value)

+                 {

+                   gmap_db_tmp_dir <- file.path(tempdir(), "gmap_db_tmp_dir")

+                   dir.create(gmap_db_tmp_dir, recursive=TRUE)

+                   cur_wd <- getwd()

+                   on.exit({unlink(gmap_db_tmp_dir, recursive=TRUE)

+                            setwd(cur_wd)})

+                   setwd(gmap_db_tmp_dir)

+

+                   db <- gmap_build(value, x, ...)

+                   

+                   db

+                 })

+

+setGeneric("snps<-", function(x, ..., value) standardGeneric("snps<-"))

+

+setReplaceMethod("snps", c("GmapGenome", "ANY"),

+                 function(x, ..., value) {

+                   snpDir <- GmapSnpDirectory(x)

+                   snps(snpDir, name = name, genome = x,

+                        iitPath = mapsDirectory(x), ...) <- value

+                   x

+                 })

+

+setGeneric("spliceSites<-",

+           function(x, ..., value) standardGeneric("spliceSites<-"))

+

+setReplaceMethod("spliceSites", c("GmapGenome", "TranscriptDb"),

+                 function(x, name, value) {

+                   exons <- exonsBy(value)

+                   exonsFlat <- unlist(exons, use.names=FALSE)

+                   exonsPart <- PartitioningByWidth(exons)

+                   exonsHead <- exonsFlat[-end(exonsPart)]

+                   donors <- flank(exonsHead, 1L, start = FALSE, both = TRUE)

+                   exonsTail <- exonsFlat[-start(exonsPart)]

+                   acceptors <- flank(exonsTail, 1L, start = TRUE, both = TRUE)

+                   sites <- c(donors, acceptors)

+                   names(sites) <- values(sites)$exon_id

+                   info <- rep(c("donor", "acceptor"), each = length(donors))

+                   info <- paste(info, unlist(width(intronsByTranscript(value)),

+                                              use.names = FALSE))

+                   values(sites) <- DataFrame(info)

+                   iit_store(sites, file.path(mapsDirectory(x), name))

+                   x

+                 })

+

+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+### Show

+###

+

+setMethod("show", "GmapGenome", function(object) {

+  cat("GmapGenome object\ngenome:", genome(object), "\ndirectory:",

+      path(directory(object)), "\n")

+})

+

+### - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -

+### Utilities

+###

+

+.normArgDb <- function(db, dir) {

+  if (!is(db, "GmapGenome"))

+    db <- GmapGenome(db, dir)

+  db

+}

+