Mercurial > repos > testtool > geo_data
view GRsetFromGEO/GRsetFromGEO.R @ 44:1edbb29139af draft
Uploaded
author | testtool |
---|---|
date | Mon, 24 Apr 2017 10:50:20 -0400 |
parents | 660897ee0d12 |
children | 7aeee7c02c4d |
line wrap: on
line source
require("GEOquery", quietly = TRUE) require("minfi", quietly = TRUE) require("FDb.InfiniumMethylation.hg19", quietly = TRUE) options(warn = -1) options("download.file.method"="wget") 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) if (is.null(GSE)) { stop("Must specify GSE") } else { options(download.file.method.GEOquery = "wget")} 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)