annotate GRsetFromGEO/GRsetFromGEO.R @ 53:46acdba0f6d8 draft

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