view GRsetFromGEO/GRsetFromGEO.R @ 35:694382fd220a draft

Uploaded

require("minfi", quietly = TRUE)

args <- commandArgs(trailingOnly = TRUE)
GSE = args[1] 
output = args[2] 

function (GSE = NULL, path = NULL, array = "IlluminaHumanMethylation450k", 
          annotation = .default.450k.annotation, what = c("Beta", "M"), 
          mergeManifest = FALSE, i = 1) 
{
  what <- match.arg(what)
  if (is.null(GSE) && is.null(path)) 
    stop("Either GSE or path must be supplied.")
  if (!is.null(GSE)) 
    gset <- GEOquery::getGEO(GSE)
  else gset <- GEOquery::getGEO(filename = file.path(path, 
                                                     list.files(path, pattern = ".soft")))
  if (length(gset) == 0) 
    stop("Empty list retrieved from GEO.")
  if (length(gset) > 1) {
    warning("More than one ExpressionSet found:\n", names(gset), 
            "\nUsing entry ", i)
    gset <- gset[[i]]
  }
  else gset <- gset[[1]]
  platform <- annotation(gset)
  if (platform != "GPL13534") 
    warning(sprintf("%s is not the platform ID associated with IlluminaHumanMethylation450k. Should be GPL13534.", 
                    platform))
  if (what == "Beta" & (min(exprs(gset)[, 1], na.rm = TRUE) < 
                        0 | max(exprs(gset)[, 1], na.rm = TRUE) > 1)) 
    warning("Values outside [0,1] detected. 'what' argument should not be Beta.")
  ann <- .getAnnotationString(c(array = array, annotation = annotation))
  if (!require(ann, character.only = TRUE)) 
    stop(sprintf("cannot load annotation package %s", ann))
  object <- get(ann)
  gr <- getLocations(object, mergeManifest = mergeManifest, 
                     orderByLocation = TRUE)
  locusNames <- names(gr)
  sampleNames(gset) <- gset$title
  common <- intersect(locusNames, fData(gset)$Name)
  if (length(common) == 0) 
    stop("No rowname matches. 'rownames' need to match IlluminaHumanMethylation450k probe names.")
  ind1 <- match(common, fData(gset)$Name)
  ind2 <- match(common, locusNames)
  preprocessing <- c(rg.norm = paste0("See GEO ", GSE, " for details"))
  if (what == "Beta") {
    out <- GenomicRatioSet(gr = gr[ind2, ], Beta = exprs(gset)[ind1, 
                                                               , drop = FALSE], M = NULL, CN = NULL, pData = pData(gset), 
                           annotation = c(array = array, annotation = annotation), 
                           preprocessMethod = preprocessing)
  }
  else {
    out <- GenomicRatioSet(gr = gr[ind2, ], Beta = NULL, 
                           M = exprs(gset)[ind1, , drop = FALSE], CN = NULL, 
                           pData = pData(gset), annotation = c(array = array, 
                                                               annotation = annotation), preprocessMethod = preprocessing)
  }
  save(out,output)
}