|
34
|
1 require("minfi", quietly = TRUE)
|
|
|
2
|
|
|
3 args <- commandArgs(trailingOnly = TRUE)
|
|
|
4 GSE = args[1]
|
|
|
5 output = args[2]
|
|
|
6
|
|
35
|
7 function (GSE = NULL, path = NULL, array = "IlluminaHumanMethylation450k",
|
|
|
8 annotation = .default.450k.annotation, what = c("Beta", "M"),
|
|
|
9 mergeManifest = FALSE, i = 1)
|
|
|
10 {
|
|
|
11 what <- match.arg(what)
|
|
|
12 if (is.null(GSE) && is.null(path))
|
|
|
13 stop("Either GSE or path must be supplied.")
|
|
|
14 if (!is.null(GSE))
|
|
|
15 gset <- GEOquery::getGEO(GSE)
|
|
|
16 else gset <- GEOquery::getGEO(filename = file.path(path,
|
|
|
17 list.files(path, pattern = ".soft")))
|
|
|
18 if (length(gset) == 0)
|
|
|
19 stop("Empty list retrieved from GEO.")
|
|
|
20 if (length(gset) > 1) {
|
|
|
21 warning("More than one ExpressionSet found:\n", names(gset),
|
|
|
22 "\nUsing entry ", i)
|
|
|
23 gset <- gset[[i]]
|
|
|
24 }
|
|
|
25 else gset <- gset[[1]]
|
|
|
26 platform <- annotation(gset)
|
|
|
27 if (platform != "GPL13534")
|
|
|
28 warning(sprintf("%s is not the platform ID associated with IlluminaHumanMethylation450k. Should be GPL13534.",
|
|
|
29 platform))
|
|
|
30 if (what == "Beta" & (min(exprs(gset)[, 1], na.rm = TRUE) <
|
|
|
31 0 | max(exprs(gset)[, 1], na.rm = TRUE) > 1))
|
|
|
32 warning("Values outside [0,1] detected. 'what' argument should not be Beta.")
|
|
|
33 ann <- .getAnnotationString(c(array = array, annotation = annotation))
|
|
|
34 if (!require(ann, character.only = TRUE))
|
|
|
35 stop(sprintf("cannot load annotation package %s", ann))
|
|
|
36 object <- get(ann)
|
|
|
37 gr <- getLocations(object, mergeManifest = mergeManifest,
|
|
|
38 orderByLocation = TRUE)
|
|
|
39 locusNames <- names(gr)
|
|
|
40 sampleNames(gset) <- gset$title
|
|
|
41 common <- intersect(locusNames, fData(gset)$Name)
|
|
|
42 if (length(common) == 0)
|
|
|
43 stop("No rowname matches. 'rownames' need to match IlluminaHumanMethylation450k probe names.")
|
|
|
44 ind1 <- match(common, fData(gset)$Name)
|
|
|
45 ind2 <- match(common, locusNames)
|
|
|
46 preprocessing <- c(rg.norm = paste0("See GEO ", GSE, " for details"))
|
|
|
47 if (what == "Beta") {
|
|
|
48 out <- GenomicRatioSet(gr = gr[ind2, ], Beta = exprs(gset)[ind1,
|
|
|
49 , drop = FALSE], M = NULL, CN = NULL, pData = pData(gset),
|
|
|
50 annotation = c(array = array, annotation = annotation),
|
|
|
51 preprocessMethod = preprocessing)
|
|
|
52 }
|
|
|
53 else {
|
|
|
54 out <- GenomicRatioSet(gr = gr[ind2, ], Beta = NULL,
|
|
|
55 M = exprs(gset)[ind1, , drop = FALSE], CN = NULL,
|
|
|
56 pData = pData(gset), annotation = c(array = array,
|
|
|
57 annotation = annotation), preprocessMethod = preprocessing)
|
|
|
58 }
|
|
|
59 save(out,output)
|
|
|
60 }
|
