# HG changeset patch # User testtool # Date 1493038259 14400 # Node ID b3761b109ca925a7e05450236ce1ad26e296fd26 # Parent 694382fd220a91dd3137d0ffd29129e92870cc62 Uploaded diff -r 694382fd220a -r b3761b109ca9 GRsetFromGEO/.Rapp.history diff -r 694382fd220a -r b3761b109ca9 GRsetFromGEO/GRsetFromGEO.R --- a/GRsetFromGEO/GRsetFromGEO.R Fri Apr 14 14:16:46 2017 -0400 +++ b/GRsetFromGEO/GRsetFromGEO.R Mon Apr 24 08:50:59 2017 -0400 @@ -1,60 +1,64 @@ 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] -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]] +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) - 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)) + + 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) + 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) + 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) } - 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) + save(out,file = output) } diff -r 694382fd220a -r b3761b109ca9 GRsetFromGEO/GRsetFromGEO.xml --- a/GRsetFromGEO/GRsetFromGEO.xml Fri Apr 14 14:16:46 2017 -0400 +++ b/GRsetFromGEO/GRsetFromGEO.xml Mon Apr 24 08:50:59 2017 -0400 @@ -26,7 +26,7 @@ - + diff -r 694382fd220a -r b3761b109ca9 GRsetFromGEO/Untitled.R --- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/GRsetFromGEO/Untitled.R Mon Apr 24 08:50:59 2017 -0400 @@ -0,0 +1,10 @@ +require("minfi", quietly = TRUE) + +args <- commandArgs(trailingOnly = TRUE) +GSE = args[1] +output = args[2] + +GRset <- getGenomicRatioSetFromGEO(GSE) + +save(GRset,output) + diff -r 694382fd220a -r b3761b109ca9 GRsetFromGEO/test-data/._.DS_Store Binary file GRsetFromGEO/test-data/._.DS_Store has changed diff -r 694382fd220a -r b3761b109ca9 GRsetFromGEO/test-data/GenomicRatioSet.Rdata Binary file GRsetFromGEO/test-data/GenomicRatioSet.Rdata has changed diff -r 694382fd220a -r b3761b109ca9 GRsetFromGEO/test-data/out.RData Binary file GRsetFromGEO/test-data/out.RData has changed