R/GMT.R

In gson: Base Class and Methods for 'gson' Format

##' parse gmt file to a data.frame
##'
##'
##' @title read.gmt
##' @rdname read-gmt
##' @param gmtfile gmt file
##' @importFrom utils stack
##' @importFrom yulab.utils yread
##' @export
##' @return data.frame
##' @author Guangchuang Yu
read.gmt <- function(gmtfile) {
  ## getGmt <- get_fun_from_pkg("GSEABase", "getGmt")
  ## geneIds <- get_fun_from_pkg("GSEABase", "geneIds")
  
  ## gmt <- getGmt(con=gmtfile)
  ## ont2gene <- geneIds(gmt) %>% stack
  ## ont2gene <- ont2gene[, c("ind", "values")]
  ## colnames(ont2gene) <- c("ont", "gene")
  
  # x <- readLines(gmtfile)
  x <- yread(gmtfile)
  ## ont2gene <- lapply(x, function(record) {
  ##     y = strsplit(record, "\t")[[1]]
  ##     data.frame(ont=y[1], gene=y[-c(1:2)])
  ## }) %>% do.call('rbind', .)
  
  
  ## first column: gene set name
  ## second column: description
  ## all the others, unequal length for genes
  res <- strsplit(x, "\t")
  names(res) <- vapply(res, function(y) y[1], character(1))
  res <- lapply(res, "[", -c(1:2))
  
  ont2gene <- stack(res)
  ont2gene <- ont2gene[, c("ind", "values")]
  colnames(ont2gene) <- c("term", "gene")
  return(ont2gene)
}

##' @rdname read-gmt
##' @param output one of 'data.frame' or 'GSON'
##' @export
read.gmt.wp <- function(gmtfile, output = "data.frame") {
  output <- match.arg(output, c("data.frame", "gson", "GSON"))
  x <- read.gmt(gmtfile)
  parts <- strsplit(as.character(x$term), "%")
  x$name <- vapply(parts, `[`, FUN.VALUE = character(1), 1)
  x$version <- vapply(parts, `[`, FUN.VALUE = character(1), 2)
  x$wpid <- vapply(parts, `[`, FUN.VALUE = character(1), 3)
  x$species <- vapply(parts, `[`, FUN.VALUE = character(1), 4)
  x$term <- NULL
  if (output == "data.frame") {
    return(x)
  }
  
  gsid2gene <- data.frame(gsid=x$wpid, gene=x$gene)
  gsid2name <- unique(data.frame(gsid=x$wpid, name=x$name))
  species <- unique(x$species)
  version <- unique(x$version)
  gson(gsid2gene = gsid2gene, 
      gsid2name = gsid2name, 
      gsname = "WikiPathways",
      species = species, 
      gsidtype = "WP",
      urlpattern = "https://www.wikipathways.org/pathways/{gsid}.html",
      version = version)
}

##' write a GSON object to GMT format
##'
##' @title write.gmt
##' @rdname read-gmt
##' @param x A `GSON` object
##' @param file output GMT file
##' @export
write.gmt <- function(x, file = "") {
  gsid2gene <- usable_memberships(x@gsid2gene)
  term2gene <- split(gsid2gene$gene, gsid2gene$gsid)
  descriptions <- setNames(names(term2gene), names(term2gene))
  if (!is.null(x@gsid2name)) {
    descriptions[x@gsid2name$gsid] <- x@gsid2name$name
  }
  res <- vapply(names(term2gene), function(gsid) {
    ## a gene set that carries no usable name is still written under its ID;
    ## the NA would otherwise reach the file as the literal string "NA"
    description <- descriptions[[gsid]]
    if (is.na(description) || !nzchar(description)) {
      description <- gsid
    }
    paste(c(gsid, description, term2gene[[gsid]]), collapse = "\t")
  }, character(1))
  names(res) <- NULL
  if (file == "") {
    return(res)
  }
  writeLines(res, file)
}