Mercurial > repos > testtool > geo_data
view GRsetFromGEO/GRsetFromGEO.R @ 36:b3761b109ca9 draft
Uploaded
author | testtool |
---|---|
date | Mon, 24 Apr 2017 08:50:59 -0400 |
parents | 694382fd220a |
children | c982fdb0e27d |
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" function (GSE = GSE, 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) }