|
52
|
1 require("BiocGenerics", quietly = TRUE)
|
|
42
|
2 require("GEOquery", quietly = TRUE)
|
|
34
|
3 require("minfi", quietly = TRUE)
|
|
36
|
4 require("FDb.InfiniumMethylation.hg19", quietly = TRUE)
|
|
34
|
5
|
|
43
|
6 options(warn = -1)
|
|
|
7 options("download.file.method"="wget")
|
|
|
8
|
|
34
|
9 args <- commandArgs(trailingOnly = TRUE)
|
|
|
10 GSE = args[1]
|
|
|
11 output = args[2]
|
|
|
12
|
|
36
|
13 getAnnotationString <- function(annotation) {
|
|
|
14 if(length(annotation) == 1)
|
|
|
15 return(sprintf("%sanno", annotation))
|
|
|
16 if(all(c("array", "annotation") %in% names(annotation)))
|
|
|
17 return(sprintf("%sanno.%s", annotation["array"], annotation["annotation"]))
|
|
|
18 stop("unable to get the annotation string for this object")
|
|
|
19 }
|
|
|
20
|
|
|
21 default.450k.annotation <- "ilmn12.hg19"
|
|
|
22
|
|
|
23
|
|
|
24
|
|
37
|
25 array = "IlluminaHumanMethylation450k"
|
|
|
26 annotation = default.450k.annotation
|
|
|
27 what = c("Beta", "M")
|
|
|
28 mergeManifest = FALSE
|
|
|
29 i = 1
|
|
|
30
|
|
43
|
31 gset <- getGEO(GSE)
|
|
|
32
|
|
|
33 if (is.null(GSE)) {
|
|
|
34 stop("Must specify GSE")
|
|
|
35 } else {
|
|
|
36 options(download.file.method.GEOquery = "wget")}
|
|
|
37
|
|
37
|
38 gset <- gset[[1]]
|
|
|
39
|
|
|
40 platform <- annotation(gset)
|
|
|
41
|
|
|
42 ann <- getAnnotationString(c(array = array, annotation = annotation))
|
|
|
43 if (!require(ann, character.only = TRUE))
|
|
|
44 stop(sprintf("cannot load annotation package %s", ann))
|
|
|
45
|
|
|
46 object <- get(ann)
|
|
|
47
|
|
|
48 gr <- getLocations(object, mergeManifest = mergeManifest,
|
|
|
49 orderByLocation = TRUE)
|
|
|
50 locusNames <- names(gr)
|
|
|
51 sampleNames(gset) <- gset$title
|
|
|
52 common <- intersect(locusNames, fData(gset)$Name)
|
|
|
53 if (length(common) == 0) {
|
|
|
54 stop("No rowname matches. 'rownames' need to match IlluminaHumanMethylation450k probe names.")
|
|
|
55 ind1 <- match(common, fData(gset)$Name)
|
|
|
56 ind2 <- match(common, locusNames)
|
|
|
57 preprocessing <- c(rg.norm = paste0("See GEO ", GSE, " for details"))
|
|
|
58 if (what == "Beta") {
|
|
|
59 out <- GenomicRatioSet(gr = gr[ind2, ], Beta = exprs(gset)[ind1,
|
|
|
60 , drop = FALSE], M = NULL, CN = NULL, pData = pData(gset),
|
|
|
61 annotation = c(array = array, annotation = annotation),
|
|
|
62 preprocessMethod = preprocessing)
|
|
35
|
63 }
|
|
37
|
64 else {
|
|
|
65 out <- GenomicRatioSet(gr = gr[ind2, ], Beta = NULL,
|
|
|
66 M = exprs(gset)[ind1, , drop = FALSE], CN = NULL,
|
|
|
67 pData = pData(gset), annotation = c(array = array,
|
|
|
68 annotation = annotation), preprocessMethod = preprocessing)
|
|
|
69 }
|
|
|
70 return(out)
|
|
35
|
71 }
|
|
53
|
72 write(out,file = output)
|
|
37
|
73
|
|
|
74
|
