R/2_geo_id.R

Defines functions get_gpl_txt get_count_txt get_ids find_anno .tinyarray_geo_download_impl geo_download .tinyarray_geo_try_annoprobe_geoChina .tinyarray_geo_manual_download_message .tinyarray_geo_read_cache_safe .tinyarray_geo_cache_is_complete .tinyarray_geo_remote_content_length .tinyarray_geo_guess_gse_from_cache_path .tinyarray_geo_annoprobe_url geo_parser

Documented in find_anno geo_download geo_parser get_count_txt get_gpl_txt get_ids

##' geo_parser
##'
##' Legacy helper for parsing GEO series matrix files.
##' Prefer [geo_download()] for normal use.
##'
##' @param gse GEO accession number
##' @param destdir The destination directory for data downloads
##' @return an ExpressionSet object list
##' @author Xiaojie Sun
##' @export
##' @examples
##' \dontrun{
##' gse = "GSE42872"
##' a = try(geo_parser(gse,destdir=tempdir()), silent = TRUE)
##' }
##' @seealso
##' \code{\link{find_anno}}

geo_parser <- function(gse,destdir = getwd()) {
  series_matrix_file = paste0(gse,'_series_matrix.txt.gz')
  if(!file.exists(paste0(destdir,"/",series_matrix_file))){
    message(paste0("Downloading ", series_matrix_file, " from GEO ..."))
  }
  if(requireNamespace("GEOquery",quietly = TRUE)) {
    eSet = tryCatch({GEOquery::getGEO(gse, destdir = destdir,getGPL = FALSE)
    },error = function(e){555})
  }else{
    eSet = 555
  }

  if(!file.exists(paste0(destdir,"/",series_matrix_file))){
    stop('Unable to download the series matrix file. Please download the series
         matrix file from GEO, place it in your working directory, and then
         rerun this function.')
  }

  if(is.list(eSet)){
    return(eSet)
  }else{
    if(!requireNamespace("Biobase",quietly = TRUE)) {
      stop("Package \"Biobase\" needed for this function to work.
           Please install it by BiocManager::install('Biobase')",call. = FALSE)
    }else{
      lines <- readLines(series_matrix_file)
      pd <- strsplit(gsub("^!|\"", "", lines), "\t", fixed = TRUE)
      width <- max(lengths(pd))
      pd <- do.call(rbind, lapply(pd, function(x) {
        length(x) <- width
        x
      }))
      k <- startsWith(pd[,1], "Sample_")
      pd <- pd[k, ]
      rownames(pd) <- sub("^Sample_", "", pd[,1])
      pd <- t(pd[,-1])
      rownames(pd) = pd[,2]
      metaData <- data.frame(labelDescription=colnames(pd))
      pd <- Biobase::AnnotatedDataFrame(data = data.frame(pd),varMetadata = metaData)
      exp <- utils::read.delim(series_matrix_file,
                               comment.char = "!",
                               row.names = 1) |>
        as.matrix()
      gpl = sub(".*(GPL\\d+).*", "\\1", lines[which(grepl("GPL", lines))[1]])
      eSet = list(Biobase::ExpressionSet(assayData = exp,phenoData = pd,annotation = gpl))
      return(eSet)
    }
  }
}

.tinyarray_geo_annoprobe_url <- function(gse) {
  gse <- toupper(as.character(gse)[1L])
  if (!grepl("^GSE[0-9]+$", gse)) {
    return(NA_character_)
  }
  down <- if (as.numeric(gsub("GSE", "", gse)) < 1000) {
    paste0("/GEOmirror/GSEnnn/", gse, "_eSet.Rdata")
  } else {
    paste0("/GEOmirror/", gsub("[0-9][0-9][0-9]$", "nnn", gse), "/", gse, "_eSet.Rdata")
  }
  paste0("http://49.235.27.111", down)
}

.tinyarray_geo_guess_gse_from_cache_path <- function(path) {
  nm <- basename(path)
  m <- regexpr("GSE[0-9]+", nm, ignore.case = TRUE)
  if (m[1L] < 0L) {
    return(NA_character_)
  }
  toupper(regmatches(nm, m)[1L])
}

.tinyarray_geo_remote_content_length <- function(url, timeout = getOption("timeout")) {
  if (!requireNamespace("httr2", quietly = TRUE)) {
    return(NA_real_)
  }
  timeout_value <- suppressWarnings(as.numeric(timeout)[1L])
  if (is.na(timeout_value) || !is.finite(timeout_value) || timeout_value <= 0) {
    timeout_value <- getOption("timeout")
  }
  timeout_value <- min(30L, max(1L, as.integer(timeout_value)))
  resp <- tryCatch(
    httr2::request(url) |>
      httr2::req_method("HEAD") |>
      httr2::req_timeout(timeout_value) |>
      httr2::req_perform(),
    error = function(e) NULL
  )
  if (is.null(resp)) {
    return(NA_real_)
  }
  cl <- httr2::resp_header(resp, "content-length")
  cl <- suppressWarnings(as.numeric(cl)[1L])
  if (is.na(cl) || !is.finite(cl) || cl <= 0) {
    return(NA_real_)
  }
  cl
}

.tinyarray_geo_cache_is_complete <- function(file,
                                             gse = NULL,
                                             cache_kind = c("rds", "rdata"),
                                             download_timeout = getOption("timeout")) {
  cache_kind <- match.arg(cache_kind)
  if (!file.exists(file)) {
    return(FALSE)
  }
  if (cache_kind == "rdata" && !is.null(gse)) {
    expected <- .tinyarray_geo_remote_content_length(
      .tinyarray_geo_annoprobe_url(gse),
      timeout = download_timeout
    )
    if (!is.na(expected)) {
      actual <- suppressWarnings(as.numeric(file.size(file))[1L])
      if (!is.na(actual) && actual != expected) {
        message(
          "The cache file is incomplete: ", file, ". ",
          "The partial file will be kept for resume."
        )
        return(FALSE)
      }
    }
  }
  TRUE
}

.tinyarray_geo_read_cache_safe <- function(file,
                                           object = NULL,
                                           gse = NULL,
                                           cache_kind = c("rds", "rdata"),
                                           download_timeout = getOption("timeout")) {
  cache_kind <- match.arg(cache_kind)
  if (!.tinyarray_geo_cache_is_complete(
    file = file,
    gse = gse,
    cache_kind = cache_kind,
    download_timeout = download_timeout
  )) {
    return(NULL)
  }
  suppressWarnings(
    tryCatch(
      .tinyarray_read_cache(file, object = object),
      error = function(e) {
        message("The cache file is incomplete: ", file, ". The partial file will be kept for resume.")
        NULL
      }
    )
  )
}

.tinyarray_geo_manual_download_message <- function(gse) {
  series_matrix_file <- paste0(gse, "_series_matrix.txt.gz")
  paste0(
    "Unable to obtain GEO data for ", gse, ". ",
    "Please download `", series_matrix_file, "` from the GEO website, ",
    "place it in your working directory, and rerun `geo_download(\"", gse, "\")`."
  )
}

.tinyarray_geo_try_annoprobe_geoChina <- function(gse,
                                                 destdir = getwd(),
                                                 download_timeout = getOption("timeout")) {
  if (!requireNamespace("curl", quietly = TRUE)) {
    return(555)
  }

  timeout_value <- suppressWarnings(as.numeric(download_timeout)[1L])
  if (is.na(timeout_value) || !is.finite(timeout_value) || timeout_value <= 0) {
    timeout_value <- getOption("timeout")
  }
  timeout_value <- max(1800L, as.integer(timeout_value))

  url <- .tinyarray_geo_annoprobe_url(gse)
  if (is.na(url) || !nzchar(url)) {
    return(555)
  }
  tpf <- file.path(destdir, paste0(gse, "_eSet.Rdata"))
  expected <- .tinyarray_geo_remote_content_length(url, timeout = timeout_value)
  last_error <- NULL

  if (!is.na(expected) && file.exists(tpf)) {
    current_size <- suppressWarnings(as.numeric(file.size(tpf))[1L])
    if (!is.na(current_size) && current_size > expected) {
      unlink(tpf)
    }
  }

  existing_size <- if (file.exists(tpf)) suppressWarnings(as.numeric(file.size(tpf))[1L]) else 0
  if (is.na(existing_size) || existing_size < 0) {
    existing_size <- 0
  }

  if (existing_size > 0) {
    message("The cache file is incomplete. Resuming download from the partial file: ", tpf)
  } else {
    message("Downloading ", basename(tpf), " from the AnnoProbe mirror ...")
  }

  handle <- curl::new_handle()
  curl::handle_setopt(
    handle,
    timeout = 0,
    connecttimeout = timeout_value,
    low_speed_time = timeout_value,
    low_speed_limit = 1,
    resume_from = existing_size
  )

  ok <- tryCatch(
    {
      curl::curl_download(
        url = url,
        destfile = tpf,
        quiet = TRUE,
        mode = if (existing_size > 0) "ab" else "wb",
        handle = handle
      )
      TRUE
    },
    error = function(e) {
      last_error <<- conditionMessage(e)
      FALSE
    }
  )

  actual_size <- if (file.exists(tpf)) suppressWarnings(as.numeric(file.size(tpf))[1L]) else 0
  if (ok) {
    if (!is.na(expected) && !is.na(actual_size) && actual_size < expected) {
      last_error <- paste0("download incomplete (expected ", expected, ", got ", actual_size, ")")
    } else {
      res <- tryCatch(
        .tinyarray_read_cache(tpf, object = c("eSet", "gset")),
        error = function(e) {
          last_error <<- conditionMessage(e)
          NULL
        }
      )
      if (!is.null(res)) {
        return(res)
      }
    }
  }

  if (file.exists(tpf) && suppressWarnings(as.numeric(file.size(tpf))[1L]) > 0) {
    if (!is.null(last_error) && grepl("404", last_error, fixed = TRUE)) {
      return(555)
    }
    warning(
      "Mirror download for ",
      gse,
      " is incomplete; the partial file was kept for resume.",
      call. = FALSE
    )
  } else if (!is.null(last_error)) {
    if (grepl("404", last_error, fixed = TRUE)) {
      return(555)
    }
    warning(
      "Mirror download failed for ",
      gse,
      ": ",
      last_error,
      call. = FALSE
    )
  }
  555
}


##' geo_download
##'
##' Download GEO expression data and return the expression matrix, phenotype
##' data, and platform information.
##'
##' Cached `eSet` objects are stored as `.rds` files in `destdir`; legacy
##' `.Rdata` cache files are still readable, and incomplete cache files will be
##' kept so the download can resume automatically.
##'
##' @param gse GEO accession number or a path to a cached `.rds` or `.RData`
##'   file containing an `ExpressionSet` or a list of `ExpressionSet` objects
##' @param by_annopbrobe whether to try the AnnoProbe cache first and fall back
##'   to GEOquery::getGEO()
##' @param simpd whether to simplify phenotype data by dropping columns with a
##'   single unique value
##' @param colon_remove whether to drop raw `characteristics_ch1*` columns
##'   while keeping the derived `:ch1` phenotype column names intact
##' @param download_timeout timeout in seconds used for GEOquery fallback and
##'   for mirror stall detection during the resumable download
##' @param destdir The destination directory for data downloads
##' @param n For data with more than one ExpressionSet, specify which one to analyze
##' @return a list with `exp`, `pd`, `gpl`, and `group_candidates`, or `NULL`
##'   when the data cannot be read or downloaded
##' @author Xiaojie Sun
##' @export
##' @examples
##' \dontrun{
##' if(requireNamespace("Biobase",quietly = TRUE)){
##'   gse = "GSE42872"
##'   a = geo_download(gse,destdir=tempdir())
##' }
##' }
##' @seealso
##' \code{\link{find_anno}}

geo_download <- function(gse, by_annopbrobe = TRUE,
                         simpd = TRUE, colon_remove = TRUE,
                         download_timeout = getOption("timeout"),
                         destdir = getwd(), n = 1) {
  tryCatch(
    .tinyarray_geo_download_impl(
      gse = gse,
      by_annopbrobe = by_annopbrobe,
      simpd = simpd,
      colon_remove = colon_remove,
      download_timeout = download_timeout,
      destdir = destdir,
      n = n
    ),
    error = function(e) {
      msg <- paste0("geo_download() failed: ", conditionMessage(e))
      tryCatch(
        warning(msg, call. = FALSE),
        error = function(w) message("Warning: ", msg)
      )
      NULL
    }
  )
}

.tinyarray_geo_download_impl <- function(gse, by_annopbrobe = TRUE,
                                         simpd = TRUE, colon_remove = TRUE,
                                         download_timeout = getOption("timeout"),
                                         destdir = getwd(), n = 1) {
  if (length(gse) != 1L || is.na(gse) || !nzchar(trimws(as.character(gse)))) {
    stop("`gse` must be one GEO accession or one cache-file path.", call. = FALSE)
  }
  gse_input <- trimws(as.character(gse))
  direct_cache <- path.expand(gse_input)
  if (!file.exists(direct_cache)) {
    cache_in_destdir <- file.path(destdir, gse_input)
    if (file.exists(cache_in_destdir)) direct_cache <- cache_in_destdir
  }

  if (file.exists(direct_cache) && !dir.exists(direct_cache)) {
    cache_kind <- if (grepl("\\.rds$", direct_cache, ignore.case = TRUE)) "rds" else "rdata"
    eSet <- .tinyarray_geo_read_cache_safe(
      direct_cache,
      object = c("eSet", "gset"),
      gse = if (cache_kind == "rdata") basename(direct_cache) else NULL,
      cache_kind = cache_kind,
      download_timeout = download_timeout
    )
    if (is.null(eSet)) {
      gse_guess <- .tinyarray_geo_guess_gse_from_cache_path(direct_cache)
      if (is.na(gse_guess)) {
        stop(
          "The cache file is incomplete and cannot be resumed because the GEO accession could not be inferred from its filename: ",
          direct_cache,
          call. = FALSE
        )
      }
      message(
        "The cache file is incomplete. The partial file will be kept and the download will resume automatically: ",
        direct_cache
      )
      return(.tinyarray_geo_download_impl(
        gse = gse_guess,
        by_annopbrobe = by_annopbrobe,
        simpd = simpd,
        colon_remove = colon_remove,
        download_timeout = download_timeout,
        destdir = dirname(direct_cache),
        n = n
      ))
    }
    gse_guess <- .tinyarray_geo_guess_gse_from_cache_path(direct_cache)
    gse <- if (!is.na(gse_guess)) gse_guess else basename(direct_cache)
  } else {
    if (grepl("\\.(rds|rdata)$", gse_input, ignore.case = TRUE)) {
      stop("Cache file does not exist: ", gse_input, call. = FALSE)
    }
    gse = toupper(gse_input)
    if (!grepl("^GSE[0-9]+$", gse)) {
      stop(
        "`gse` must be a GEO accession such as \"GSE7305\" or an existing cache-file path.",
        call. = FALSE
      )
    }
    if (by_annopbrobe) {
      cache_file = file.path(destdir,paste0(gse,"_eSet.rds"))
      legacy_cache_file = file.path(destdir,paste0(gse,"_eSet.Rdata"))
      eSet <- NULL
      if (file.exists(cache_file)) {
        eSet <- .tinyarray_geo_read_cache_safe(
          cache_file,
          object = "eSet",
          cache_kind = "rds",
          download_timeout = download_timeout
        )
      }
      if (is.null(eSet) && file.exists(legacy_cache_file)) {
        eSet <- .tinyarray_geo_read_cache_safe(
          legacy_cache_file,
          object = c("eSet", "gset"),
          gse = gse,
          cache_kind = "rdata",
          download_timeout = download_timeout
        )
        if (!is.null(eSet)) {
          saveRDS(eSet, cache_file)
          if (file.exists(legacy_cache_file)) unlink(legacy_cache_file)
        }
      }
      if (is.null(eSet)) {
        eSet <- .tinyarray_geo_try_annoprobe_geoChina(
          gse = gse,
          destdir = destdir,
          download_timeout = download_timeout
        )

        if (!is.list(eSet)) {
          if (file.exists(legacy_cache_file) && file.size(legacy_cache_file) > 0) {
            stop(
              "The AnnoProbe mirror download for ",
              gse,
              " did not finish. The partial cache was kept for resume, so rerun ",
              "`geo_download()` to continue the same download.",
              call. = FALSE
            )
          }
          message("This data is not indexed by AnnoProbe, downloading from GEO ...")
          if (!requireNamespace("GEOquery", quietly = TRUE)) {
            stop("Package \"GEOquery\" needed for this function to work.
         Please install it by BiocManager::install('GEOquery')", call. = FALSE)
          }
          timeout_value <- suppressWarnings(as.numeric(download_timeout)[1L])
          if (is.na(timeout_value) || !is.finite(timeout_value) || timeout_value <= 0) {
            timeout_value <- getOption("timeout")
          }
          timeout_value <- max(1, as.integer(timeout_value))
          old_timeout <- options(timeout = timeout_value)
          eSet <- tryCatch({GEOquery::getGEO(gse, destdir = destdir, getGPL = FALSE)
          }, error = function(e) {555})
          options(old_timeout)
          if (!is.list(eSet)) {
            stop(.tinyarray_geo_manual_download_message(gse), call. = FALSE)
          }
        }
        if (is.list(eSet)) {
          saveRDS(eSet, cache_file)
          if (file.exists(legacy_cache_file)) unlink(legacy_cache_file)
        }
      }
    } else {
      if(!requireNamespace("GEOquery",quietly = TRUE)) {
        stop("Package \"GEOquery\" needed for this function to work.
         Please install it by BiocManager::install('GEOquery')",call. = FALSE)
      }
      timeout_value <- suppressWarnings(as.numeric(download_timeout)[1L])
    if (is.na(timeout_value) || !is.finite(timeout_value) || timeout_value <= 0) {
      timeout_value <- getOption("timeout")
    }
      timeout_value <- max(1, as.integer(timeout_value))
      old_timeout <- options(timeout = timeout_value)
      eSet <- tryCatch({GEOquery::getGEO(gse, destdir = destdir,getGPL = FALSE)
      },error = function(e){555})
      options(old_timeout)
      if (!is.list(eSet)) {
        stop(.tinyarray_geo_manual_download_message(gse), call. = FALSE)
      }
    }
  }
  if (inherits(eSet, "ExpressionSet")) eSet <- list(eSet)
  if(!is.list(eSet)) {
    stop(.tinyarray_geo_manual_download_message(gse), call. = FALSE)
  }
  if(length(n)!=1) stop("only one ExpresssionSet can be analyzed")
  if(length(eSet)==1 & n!=1) {
    n = 1
    warning("this data only has one ExpressionSet object")
  }
  if(!requireNamespace("Biobase",quietly = TRUE)) {
    stop("Package \"Biobase\" needed for this function to work.
         Please install it by BiocManager::install('Biobase')",call. = FALSE)
  }
  exp <- Biobase::exprs(eSet[[n]])
  pd <- Biobase::pData(eSet[[n]])
  if(simpd){
    colname <- vector("character")
    count <- vector("integer")
    for (i in 1:ncol(pd)) {
      colname[i] = colnames(pd)[[i]]
      count[i] = nrow(pd[!duplicated(pd[, i]), ])
    }
      df <- data.frame(colname, count, stringsAsFactors = FALSE)
      df <- df[order(df$count, decreasing = TRUE), , drop = FALSE]
      df <- df[df$count > 1, , drop = FALSE]
      pd <-  pd[, df$colname, drop = FALSE]
      pd <- pd[, !(colnames(pd) %in% c("geo_accession", "supplementary_file")), drop = FALSE]
  }
  pd <- .tinyarray_geo_drop_redundant_columns(
    pd,
    drop_url = TRUE,
    drop_characteristics = colon_remove
  )
  drop_pd_cols <- intersect(colnames(pd), c("last_update_date", "relation", "relation."))
  if (length(drop_pd_cols)) {
    pd <- pd[, setdiff(colnames(pd), drop_pd_cols), drop = FALSE]
  }
  if (ncol(pd) == 0L) {
    pd_full <- Biobase::pData(eSet[[n]])
    pd_full <- .tinyarray_geo_drop_redundant_columns(
      pd_full,
      drop_url = TRUE,
      drop_characteristics = colon_remove
    )
    drop_pd_full_cols <- intersect(colnames(pd_full), c("last_update_date", "relation", "relation."))
    pd_full <- pd_full[, setdiff(colnames(pd_full), drop_pd_full_cols), drop = FALSE]
    if (ncol(pd_full) == 0L) {
      pd <- pd_full[, 0, drop = FALSE]
      warning("pd simplification removed all usable columns; returning an empty pd.", call. = FALSE)
    } else {
      pd <- pd_full[, 1, drop = FALSE]
      warning("pd simplification removed all columns; keeping the first usable column.", call. = FALSE)
    }
  }
  p1 = identical(rownames(pd),colnames(exp))
  p2 = all(rownames(pd) %in% colnames(exp) & colnames(exp) %in% rownames(pd))
  if(!p1) {
    exp = exp[,match(rownames(pd),colnames(exp))]
    if(!p2) {
      exp = exp[,intersect(rownames(pd),colnames(exp))]
      pd = pd[intersect(rownames(pd),colnames(exp)),]
    }
  }
  gpl <- eSet[[n]]@annotation
  pd2 <- as.data.frame(lapply(pd, as.character), stringsAsFactors = FALSE, check.names = FALSE)
  drop_pd2_cols <- intersect(colnames(pd2), c("last_update_date", "relation", "relation."))
  if (length(drop_pd2_cols)) {
    pd2 <- pd2[, setdiff(colnames(pd2), drop_pd2_cols), drop = FALSE]
  }
  rownames(pd2) <- rownames(pd)
  group_candidates <- .tinyarray_geo_group_candidates(
    pd2,
    group_col = "title",
    search_other_cols = TRUE
  )
  re = list(exp=exp,pd=pd2,gpl=gpl,group_candidates=group_candidates)
  if(!is.null(dim(exp)) && nrow(exp)!=0 && any(is.na(exp)|is.nan(exp))) {
    warning("NA or NaN values detected")
  }else if (any(exp<0)) {
    warning("negative values detected")
  } else if(!is.null(dim(exp)) && nrow(exp)!=0){
    message(paste(nrow(exp),"probes,",
                  ncol(exp),"samples",
                  "from",min(exp),
                  "to",max(exp)))}
  return(re)
}

##' Find an annotation table for a GEO platform accession.
##'
##' The function prints a short explanation plus the code used to obtain
##' annotations. It only tests AnnoProbe `type` values when no Bioconductor
##' annotation package source exists for the GPL accession.
##'
##' @param gpl a GEO platform accession
##' @param install whether to install and library the package
##' @param update whether to update the package
##' @return No return value. The function prints available annotation resources
##'   and optionally installs a Bioconductor annotation package.
##' @author Xiaojie Sun
##' @export
##' @examples
##' \dontrun{
##' try(find_anno("GPL570"), silent = TRUE)
##' }
##' @seealso
##' \code{\link{geo_download}}

find_anno <-function(gpl,install = FALSE,update = FALSE){
  gpl = toupper(gpl)
  bioc_msg <- NULL
  anno_msg <- NULL

  if(gpl %in% pkg_all$gpl) {
    qz = pkg_all$bioc_package[pkg_all$gpl== gpl]
    bioc_msg <- .tinyarray_bioc_annotation_code(qz)
    if(install){
      if(!requireNamespace("BiocManager", quietly = TRUE)) {
        stop("Package \"BiocManager\" needed for this function to work. Please install it by install.packages('BiocManager')",
             call. = FALSE)
      }
      if(!suppressMessages(requireNamespace(paste0(qz,".db")))){
        BiocManager::install(paste0(qz,".db"),update = update,ask = FALSE)
        suppressMessages(requireNamespace(paste0(qz,".db")))
      }
    }
  }

  if (is.null(bioc_msg) && requireNamespace("AnnoProbe", quietly = TRUE)) {
    anno_try <- .tinyarray_try_annoprobe_idmap(gpl)
    if (!is.null(anno_try) && !is.null(anno_try$code)) {
      anno_msg <- anno_try$code
    }
  }

  if (!is.null(bioc_msg)) {
    message("The following code can be used to obtain annotations:")
    message(bioc_msg)
  }
  if (!is.null(anno_msg)) {
    if (!is.null(bioc_msg)) {
      message("The following code snippet is also available:")
    } else {
      message("The following code can be used to obtain annotations:")
    }
    message(anno_msg)
  }
  if (is.null(bioc_msg) && is.null(anno_msg)) {
    message(
      "No annotation table is available for ",
      gpl,
      ". Please extract it manually from the GPL platform table. ",
      "You can use `get_gpl_txt(\"",
      gpl,
      "\")` to get the download link for the GPL platform table."
    )
  }
}

##' Get probe-to-symbol ids for a GEO platform accession.
##'
##' Prefer the available Bioconductor annotation package when possible, and
##' fall back to `AnnoProbe::idmap()` when no Bioconductor annotation package
##' source exists. The function prints a short explanation plus the code it
##' used.
##'
##' @param gpl a GEO platform accession
##' @param type optional AnnoProbe `type` setting. If omitted, the function
##'   tries candidate types automatically.
##' @param install whether to install and library the package
##' @param update whether to update the package
##' @return a cleaned two-column data.frame with `probe_id` and `symbol`
##' @author Xiaojie Sun
##' @export
##' @examples
##' \dontrun{
##' get_ids("GPL570")
##' }
##' @seealso
##' \code{\link{find_anno}}

get_ids <- function(gpl, type = NULL, install = FALSE, update = FALSE) {
  gpl <- toupper(gpl)
  anno <- NULL
  anno_type <- NULL
  anno_code <- NULL
  has_bioc_source <- gpl %in% pkg_all$gpl
  bioc_label <- NULL

  if (has_bioc_source) {
    qz <- pkg_all$bioc_package[pkg_all$gpl == gpl]
    bioc_label <- paste0(qz, ".db")
    if (install) {
      if (!requireNamespace("BiocManager", quietly = TRUE)) {
        stop(
          "Package \"BiocManager\" needed for this function to work. Please install it by install.packages('BiocManager')",
          call. = FALSE
        )
      }
      if (!suppressMessages(requireNamespace(bioc_label, quietly = TRUE))) {
        BiocManager::install(bioc_label, update = update, ask = FALSE)
      }
    }
    if (requireNamespace(bioc_label, quietly = TRUE)) {
      anno <- tryCatch(
        AnnotationDbi::toTable(
          getExportedValue(bioc_label, paste0(qz, "SYMBOL"))
        ),
        error = function(e) {
          warning(
            "Failed to read ",
            bioc_label,
            ": ",
            conditionMessage(e),
            call. = FALSE
          )
          NULL
        }
      )
      if (!is.null(anno)) {
        anno_type <- "bioc"
        anno_code <- .tinyarray_bioc_annotation_code(qz)
      }
    } else {
      warning(
        "A Bioconductor annotation package exists for ",
        gpl,
        " but ",
        bioc_label,
        " is not installed; falling back to AnnoProbe.",
        call. = FALSE
      )
    }
  }

  if (!is.null(anno) && !is.null(type)) {
    warning(
      "`type` is ignored because a Bioconductor annotation package source exists.",
      call. = FALSE
    )
  }

  if (is.null(anno)) {
    if (!requireNamespace("AnnoProbe", quietly = TRUE)) {
      warning(
        "Package \"AnnoProbe\" is not available, so the fallback cannot be used for ",
        gpl,
        ".",
        call. = FALSE
      )
      return(NULL)
    }
    requested_type <- if (!is.null(type) && length(type)) unique(as.character(type)) else NULL
    anno_try <- .tinyarray_try_annoprobe_idmap(gpl, type = requested_type, destdir = tempdir())
    if (!is.null(anno_try) && !is.null(anno_try$data)) {
      anno <- anno_try$data
      anno_type <- anno_try$type
      anno_code <- anno_try$code
      if (!is.null(requested_type) && length(requested_type) && !identical(requested_type[[1L]], anno_type)) {
        warning(
          "AnnoProbe::idmap(type = \"",
          requested_type[[1L]],
          "\") failed for ",
          gpl,
          "; using type = \"",
          anno_type,
          "\" instead.",
          call. = FALSE
        )
      }
    } else {
      warning(
        "AnnoProbe::idmap() failed for ",
        gpl,
        " with the available `type` settings. Please use `get_gpl_txt(\"",
        gpl,
        "\")` to get the download link for the GPL platform table.",
        call. = FALSE
      )
      return(NULL)
    }
  }
  anno <- .tinyarray_standardize_probe_symbol_table(anno)
  if (is.null(anno) || !nrow(anno)) {
    message(
      "No annotation table is available for ",
      gpl,
      ". Please extract it manually from the GPL platform table. ",
      "You can use `get_gpl_txt(\"",
      gpl,
      "\")` to get the download link for the GPL platform table."
    )
    return(NULL)
  }
  if (!is.null(anno_code)) {
    message("The following code was used to obtain the annotation table:")
    message(anno_code)
  }
  attr(anno, "source_type") <- anno_type
  anno
}

##' get count from GEO
##'
##' get RNA-seq count file from GEO database
##'
##' @inheritParams geo_download
##' @param download download the txt file or not
##' @return No return value. When `download = TRUE`, the count file is saved in
##'   `destdir`.
##' @author Xiaojie Sun
##' @export
##' @examples
##' get_count_txt("GSE162550",destdir = tempdir())
##' @seealso
##' \code{\link{geo_download}}

get_count_txt <- function(gse,destdir = getwd(),download = FALSE){
  if(!grepl("^(GSE|gse)\\d+$", gse))stop("wrong GSE accession")
  gse = toupper(gse)
  url = paste0("https://www.ncbi.nlm.nih.gov/geo/download/?type=rnaseq_counts&acc=",
             gse,
             "&format=file&file=",
             gse,
             "_raw_counts_GRCh38.p13_NCBI.tsv.gz")
  message(url)

  if(download){utils::download.file(url,destfile = file.path(destdir,paste0(gse,"_raw_counts_GRCh38.p13_NCBI.tsv.gz")))
    message("If the download fails, check that your data is RNA-seq data.")}
}


##' get gpl txt from GEO
##'
##' get gpl annotation txt file from GEO database
##'
##' @inheritParams geo_download
##' @inheritParams get_count_txt
##' @param gpl gpl accession from GEO database
##' @return No return value. When `download = TRUE`, the GPL annotation file is
##'   saved in `destdir`.
##' @author Xiaojie Sun
##' @export
##' @examples
##' get_gpl_txt("GPL23270",destdir = tempdir())
##' @seealso
##' \code{\link{geo_download}}

get_gpl_txt = function(gpl,destdir = getwd(),download = FALSE){
  if(!grepl("^(GPL|gpl)\\d+$", gpl))stop("wrong GPL accession")
  gpl = toupper(gpl)
  url = paste0("https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=",
               gpl,
               "&targ=self&form=text&view=data")
  message(url)
  if(download){utils::download.file(url,destfile = file.path(destdir,paste0(gpl,".txt")))
    message("If the download fails, check that your data is microarray data.")}
}


utils::globalVariables(c("pkg_all","exists_anno_list","gset"))

Try the tinyarray package in your browser

Any scripts or data that you put into this service are public.

tinyarray documentation built on Aug. 2, 2026, 9:07 a.m.