Mercurial > repos > testtool > geo_data

comparison GRsetFromGEO/GRsetFromGEO.R @ 36:b3761b109ca9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author testtool
date Mon, 24 Apr 2017 08:50:59 -0400
parents 694382fd220a
children c982fdb0e27d
comparison
equal deleted inserted replaced
35:694382fd220a 36:b3761b109ca9
1 require("minfi", quietly = TRUE) 1 require("minfi", quietly = TRUE)
2 require("BiocGenerics", quietly = TRUE)
3 require("data.table", quietly = TRUE)
4 require("GEOquery", quietly = TRUE)
5 require("rtracklayer", quietly = TRUE)
6 require("FDb.InfiniumMethylation.hg19", quietly = TRUE)
2 7
3 args <- commandArgs(trailingOnly = TRUE) 8 args <- commandArgs(trailingOnly = TRUE)
4 GSE = args[1] 9 GSE = args[1]
5 output = args[2] 10 output = args[2]
6 11
7 function (GSE = NULL, path = NULL, array = "IlluminaHumanMethylation450k", 12 getAnnotationString <- function(annotation) {
8 annotation = .default.450k.annotation, what = c("Beta", "M"), 13 if(length(annotation) == 1)
9 mergeManifest = FALSE, i = 1) 14 return(sprintf("%sanno", annotation))
10 { 15 if(all(c("array", "annotation") %in% names(annotation)))
11 what <- match.arg(what) 16 return(sprintf("%sanno.%s", annotation["array"], annotation["annotation"]))
12 if (is.null(GSE) && is.null(path)) 17 stop("unable to get the annotation string for this object")
13 stop("Either GSE or path must be supplied.") 18 }
14 if (!is.null(GSE)) 19
15 gset <- GEOquery::getGEO(GSE) 20 default.450k.annotation <- "ilmn12.hg19"
16 else gset <- GEOquery::getGEO(filename = file.path(path, 21
17 list.files(path, pattern = ".soft"))) 22
18 if (length(gset) == 0) 23
19 stop("Empty list retrieved from GEO.") 24 function (GSE = GSE, array = "IlluminaHumanMethylation450k", annotation = default.450k.annotation,
20 if (length(gset) > 1) { 25 what = c("Beta", "M"), mergeManifest = FALSE, i = 1){
21 warning("More than one ExpressionSet found:\n", names(gset), 26
22 "\nUsing entry ", i) 27
23 gset <- gset[[i]] 28 gset <- getGEO(GSE)
24 } 29 gset <- gset[[1]]
25 else gset <- gset[[1]] 30
26 platform <- annotation(gset) 31 platform <- annotation(gset)
27 if (platform != "GPL13534") 32
28 warning(sprintf("%s is not the platform ID associated with IlluminaHumanMethylation450k. Should be GPL13534.", 33 ann <- getAnnotationString(c(array = array, annotation = annotation))
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)) 34 if (!require(ann, character.only = TRUE))
35 stop(sprintf("cannot load annotation package %s", ann)) 35 stop(sprintf("cannot load annotation package %s", ann))
36
36 object <- get(ann) 37 object <- get(ann)
38
37 gr <- getLocations(object, mergeManifest = mergeManifest, 39 gr <- getLocations(object, mergeManifest = mergeManifest,
38 orderByLocation = TRUE) 40 orderByLocation = TRUE)
39 locusNames <- names(gr) 41 locusNames <- names(gr)
40 sampleNames(gset) <- gset$title 42 sampleNames(gset) <- gset$title
41 common <- intersect(locusNames, fData(gset)$Name) 43 common <- intersect(locusNames, fData(gset)$Name)
42 if (length(common) == 0) 44 if (length(common) == 0) {
43 stop("No rowname matches. 'rownames' need to match IlluminaHumanMethylation450k probe names.") 45 stop("No rowname matches. 'rownames' need to match IlluminaHumanMethylation450k probe names.")
44 ind1 <- match(common, fData(gset)$Name) 46 ind1 <- match(common, fData(gset)$Name)
45 ind2 <- match(common, locusNames) 47 ind2 <- match(common, locusNames)
46 preprocessing <- c(rg.norm = paste0("See GEO ", GSE, " for details")) 48 preprocessing <- c(rg.norm = paste0("See GEO ", GSE, " for details"))
47 if (what == "Beta") { 49 if (what == "Beta") {
48 out <- GenomicRatioSet(gr = gr[ind2, ], Beta = exprs(gset)[ind1, 50 out <- GenomicRatioSet(gr = gr[ind2, ], Beta = exprs(gset)[ind1,
49 , drop = FALSE], M = NULL, CN = NULL, pData = pData(gset), 51 , drop = FALSE], M = NULL, CN = NULL, pData = pData(gset),
50 annotation = c(array = array, annotation = annotation), 52 annotation = c(array = array, annotation = annotation),
51 preprocessMethod = preprocessing) 53 preprocessMethod = preprocessing)
54 }
55 else {
56 out <- GenomicRatioSet(gr = gr[ind2, ], Beta = NULL,
57 M = exprs(gset)[ind1, , drop = FALSE], CN = NULL,
58 pData = pData(gset), annotation = c(array = array,
59 annotation = annotation), preprocessMethod = preprocessing)
60 }
61 return(out)
52 } 62 }
53 else { 63 save(out,file = output)
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 } 64 }