view GRsetFromGEO/GRsetFromGEO.R @ 39:925945cca45f draft

Uploaded
author testtool
date Mon, 24 Apr 2017 09:29:38 -0400
parents c982fdb0e27d
children 234e990e8e1d
line wrap: on
line source

require("minfi", quietly = TRUE)
require("BiocGenerics", quietly = TRUE)
require("data.table", quietly = TRUE)
require("GEOquery", quietly = TRUE)
require("rtracklayer", quietly = TRUE)
require("FDb.InfiniumMethylation.hg19", quietly = TRUE)

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

getAnnotationString <- function(annotation) {
  if(length(annotation) == 1)
    return(sprintf("%sanno", annotation))
  if(all(c("array", "annotation") %in% names(annotation)))
    return(sprintf("%sanno.%s", annotation["array"], annotation["annotation"]))
  stop("unable to get the annotation string for this object")
}

default.450k.annotation <- "ilmn12.hg19"



array = "IlluminaHumanMethylation450k"
annotation = default.450k.annotation
what = c("Beta", "M")
mergeManifest = FALSE
i = 1

gset <- getGEO(GSE)
gset <- gset[[1]]

platform <- annotation(gset)

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)
  }
  return(out)
}
save(out,file = output)