Mercurial > repos > hackdna > pathprint

comparison galaxy-pathprint.r @ 0:0cebe436a553 draft default tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Initial upload.
author hackdna
date Fri, 17 May 2013 14:29:33 -0400
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:0cebe436a553
1
2 ## List of arrays supported by PathPrint as of version 1.2.2 This will need to
3 ## be updated/replaced with a better version as PathPrint's support changes.
4 supportedArrays <- function() {
5 return(matrix(c("GPL72", "DrosGenome1", "Affymetrix Drosophila Genome Array", "drosophila",
6 "GPL85", "RG_U34A", "Affymetrix Rat Genome U34 Array", "rat",
7 "GPL91", "HG_U95A", "Affymetrix Human Genome U95A Array", "human",
8 "GPL96", "HG-U133A", "Affymetrix Human Genome U133A Array", "human",
9 "GPL200", "Celegans", "Affymetrix C. elegans Genome Array", "C.elegans",
10 "GPL339", "MOE430A", "Affymetrix Mouse Expression 430A Array", "mouse",
11 "GPL341", "RAE230A", "Affymetrix Rat Expression 230A Array", "rat",
12 "GPL570", "HG-U133_Plus_2", "Affymetrix Human Genome U133 Plus 2.0 Array", "human",
13 "GPL571", "HG-U133A_2", "Affymetrix Human Genome U133A 2.0 Array", "human",
14 "GPL1261", "Mouse430_2", "Affymetrix Mouse Genome 430 2.0 Array", "mouse",
15 "GPL1319", "Zebrafish", "Affymetrix Zebrafish Genome Array", "zebrafish",
16 "GPL1322", "Drosophila_2", "Affymetrix Drosophila Genome 2.0 Array", "drosophila",
17 "GPL1355", "Rat230_2", "Affymetrix Rat Genome 230 2.0 Array", "rat",
18 "GPL2700", "", "Sentrix HumanRef-8 Expression BeadChip", "human",
19 "GPL2986", "", "ABI Human Genome Survey Microarray Version 2", "human",
20 "GPL2995", "", "ABI Mouse Genome Survey Microarray", "mouse",
21 "GPL3921", "HT_HG-U133A", "Affymetrix HT Human Genome U133A Array", "human",
22 "GPL4685", "U133AAofAv2", "Affymetrix GeneChip HT-HG_U133A Early Access Array", "human",
23 "GPL6102", "", "Illumina human-6 v2.0 expression beadchip", "human",
24 "GPL6103", "", "Illumina mouseRef-8 v1.1 expression beadchip", "mouse",
25 "GPL6104", "", "Illumina humanRef-8 v2.0 expression beadchip", "human",
26 "GPL6105", "", "Illumina mouse-6 v1.1 expression beadchip", "mouse",
27 "GPL6333", "", "Illumina Mouse Ref-6 V1", "mouse",
28 "GPL6883", "", "Illumina HumanRef-8 v3.0 expression beadchip", "human",
29 "GPL6884", "", "Illumina HumanWG-6 v3.0 expression beadchip", "human",
30 "GPL6885", "", "Illumina MouseRef-8 v2.0 expression beadchip", "mouse",
31 "GPL6887", "", "Illumina MouseWG-6 v2.0 expression beadchip", "mouse",
32 "GPL6947", "", "Illumina HumanHT-12 V3.0 expression beadchip", "human",
33 "GPL8300", "HG_U95Av2", "Affymetrix Human Genome U95 Version 2 Array", "human",
34 "GPL8321", "Mouse430A_2", "Affymetrix Mouse Genome 430A 2.0 Array", "mouse"),
35 ncol = 4, byrow = TRUE))
36 }
37
38 fingerprintGSM <- function(geoID) {
39 if (is.existingGeoID(geoID)) {
40 ## If we've already got the GEO ID in the matrix, don't bother loading it
41 ## again, just pull out the data we need.
42 gsm.fingerprint <- GEO.fingerprint.matrix[,geoID]
43 } else {
44 ## We weren't able to find the GEO ID in the matrix so try to load all the
45 ## data we need from NCBI.
46 try(gsm <- getGEO(geoID, GSElimits = NULL), silent = FALSE)
47 if (! exists("gsm")) {
48 sink(stderr(), type = "message");
49 stop(sprintf("Unable to load GEO id '%s': is the GEO id valid?\n", geoID))
50 }
51
52 gsm.exprs <- Table(gsm)
53 gsm.platform <- Meta(gsm)$platform_id
54 gsm.species <- Meta(gsm)$organism_ch1
55
56 ## Compute the fingerprint for the GSM id.
57 gsm.fingerprint <- exprs2fingerprint(exprs = gsm.exprs,
58 platform = gsm.platform,
59 species = gsm.species,
60 progressBar = FALSE)
61 }
62
63 ## Convert the fingerprint to a data frame.
64 gsm.fingerprint <- as.data.frame(gsm.fingerprint)
65
66 ## Tag the fingerprint as having been calculated on a GSM id.
67 attr(gsm.fingerprint, 'fingerprintType') <- 'GSM'
68 attr(gsm.fingerprint, 'fingerprintGEO') <- geoID
69
70 return(gsm.fingerprint)
71 }
72
73 fingerprintGSE <- function(geoID) {
74 if (is.existingGeoID(geoID) && 0) {
75 ## If we've already got the GEO ID in the matrix, don't bother loading it
76 ## again, just pull out the data we need.
77 gse <- na.omit(GEO.metadata.matrix[GEO.metadata.matrix$GSE == geoID,])
78 gse.gsm <- gse[,"GSM"]
79
80 gse.fingerprint <- GEO.fingerprint.matrix[,gse.gsm]
81 } else {
82 try(gse <- getGEO(geoID), silent = FALSE)
83 if (! exists("gse")) {
84 sink(stderr(), type = "message");
85 stop(sprintf("Unable to load GEO id '%s': is the GEO id valid?\n", geoID))
86 } else {
87 print(gse)
88 }
89
90 gse.exprs <- exprs(gse[[1]])
91 gse.platform <- annotation(gse[[1]])
92 gse.species <- as.character(unique(phenoData(gse[[1]])$organism_ch1))
93
94 ## Compute the fingerprint for the GSE id.
95 gse.fingerprint <- exprs2fingerprint(exprs = gse.exprs,
96 platform = gse.platform,
97 species = gse.species,
98 progressBar = FALSE)
99 }
100
101 ## Tag the fingerprint as having been calculated on a GSE id.
102 attr(gse.fingerprint, 'fingerprintType') <- 'GSE'
103 attr(gse.fingerprint, 'fingerprintGEO') <- geoID
104
105 return(gse.fingerprint)
106 }
107
108 generateFingerprint <- function(geoID) {
109 if (is.geoID(geoID)) {
110 if (is.geoGSM(geoID)) {
111 return(fingerprintGSM(geoID))
112 } else {
113 return(fingerprintGSE(geoID))
114 }
115 } else {
116 sink(stderr(), type = "message");
117 stop(paste("not a GEO ID:", geoID))
118 }
119 }
120
121 loadFingerprint <- function(file) {
122 fingerprint <- read.delim(file, header = TRUE, sep = "\t");
123 return(fingerprint);
124 }
125
126 saveFingerprint <- function(fingerprint, file) {
127 if (! is.null(attr(fingerprint, 'fingerprintType'))) {
128 if ((attr(fingerprint, 'fingerprintType') == 'GSE')) {
129 data <- fingerprint[,1:ncol(fingerprint)]
130 cols <- colnames(fingerprint);
131
132 write.table(data, row.names = TRUE,
133 col.names = cols,
134 file = file,
135 quote = FALSE,
136 sep = "\t")
137 } else {
138 data <- data.frame(fingerprint[,1])
139 rownames(data) <- rownames(fingerprint)
140 write.table(data, row.names = TRUE,
141 col.names = attr(fingerprint, 'fingerprintGEO'),
142 file = file,
143 quote = FALSE,
144 sep = "\t")
145 }
146 } else {
147 sink(stderr(), type = "message");
148 stop("Unable to save fingerprint: unknown fingerprint type")
149 }
150 }
151
152 saveConsensus <- function(consensus, file) {
153 write.table(consensus, row.names = TRUE,
154 col.names = TRUE,
155 file = file,
156 quote = FALSE,
157 sep = "\t")
158 }
159
160 calculateDistanceToGEO <- function(consensus) {
161 # sample from matrix to speed it up
162 # sample <- sample(dim(GEO.fingerprint.matrix)[2], 10000, replace=FALSE)
163 # GEO.distance <- consensusDistance(consensus, GEO.fingerprint.matrix[,sample])
164 GEO.distance <- consensusDistance(consensus, GEO.fingerprint.matrix)
165
166 similar.GEO <- GEO.metadata.matrix[
167 match(rownames(GEO.distance), GEO.metadata.matrix$GSM),
168 c("GSM", "GSE", "GPL", "Source")]
169
170 similar.GEO <- cbind(similar.GEO[1:sum(GEO.distance$pvalue < 0.01),],
171 GEO.distance[1:sum(GEO.distance$pvalue < 0.01),])
172 return(similar.GEO)
173
174 }
175
176 calculateDistanceToPluripotent <- function(fingerprint, verbose = 0) {
177 #if (ncol(fingerprint) > 1) {
178 if (verbose) print("Calculating distances from pluripotent consensus for multiple sample")
179 ## Calculate the distance from the pluripotent consensus for a multiple sample.
180 pluripotent.consensus <- consensusFingerprint(GEO.fingerprint.matrix[,pluripotents.frame$GSM], 0.9)
181 fingerprint.distance <- consensusDistance(pluripotent.consensus, fingerprint)
182 if(ncol(fingerprint)==1) {
183 rownames(fingerprint.distance) <- colnames(fingerprint)
184 }
185 #} else {
186 # if (verbose) print("Calculating distances from pluripotent consensus for single sample")
187 # ## Calculate the distance from the pluripotent consensus for a single sample.
188 # pluripotent <- GEO.fingerprint.matrix[,pluripotents.frame$GSM]
189 # fingerprint.distance <- data.frame(consensusDistance(fingerprint, pluripotent))
190 #}
191
192 return(fingerprint.distance)
193 }
194
195 saveDistance <- function(distance, file, pvalue = 1) {
196 write.table(distance, row.names = TRUE,
197 col.names = TRUE,
198 file = file,
199 quote = FALSE,
200 sep = "\t")
201 }
202
203 loadDistance <- function(file, pvalue = 1) {
204 distance <- read.delim(file, header = TRUE, sep = "\t");
205 return(distance);
206 }
207
208 generateHistograms <- function(distance, filename) {
209 if (nrow(distance) == 0) {
210 sink(stderr(), type = "message");
211 stop("Unable to generate histogram: no distances to plot - try increasing p-value?")
212 }
213
214 pdf(filename, width = 8, height = 8)
215
216 xlab1 <- "Distance of GEO records from pluripotent consensus"
217 xlab2 <- "Distance of submitted data from pluripotent consensus"
218
219 ## Calculate the distance of the pluripotent to GEO.
220 pluripotent.consensus <- consensusFingerprint(GEO.fingerprint.matrix[,pluripotents.frame$GSM], 0.9)
221 geo.distance = calculateDistanceToGEO(pluripotent.consensus)
222
223 ## Output the first histogram.
224 par(oma = c(5, 2, 2, 2), mfcol = c(2,1), mar = c(0, 4, 4, 2))
225 hist(geo.distance[,"distance"], col ="grey",
226 main = "",
227 nclass = 50,
228 xlab = xlab1,
229 xlim = c(0,1))
230
231 ## Output the second histogram.
232 par(mar = c(7, 4, 4, 2))
233 hist(distance[,"distance"], cex.lab = 0.8,
234 col = "green",
235 main = "",
236 xlab = xlab2,
237 xlim = c(0,1))
238 mtext(cex = 0.8,
239 line = 0.5,
240 side = 3,
241 text = "Distance of GEO records from pluripotent consensus")
242
243 invisible(dev.off())
244 }
245
246 generateHeatmap <- function(fingerprint, sdev, filename) {
247
248 if(ncol(fingerprint) <= 1) {
249 sink(stderr(), type = "message");
250 stop("unable to generate heatmap on fingerprint with single column. Requires a dataset with multiple samples.",
251 call. = FALSE)
252 }
253
254 heatmap.data <- fingerprint[apply(fingerprint, 1, sd) > sdev, ]
255 if(dim(heatmap.data)[1] < 2) {
256 sink(stderr(), type="message");
257 stop("unable to generate heatmap - try lowering the standard deviation cutoff",
258 call. = FALSE)
259 }
260
261 library(pheatmap)
262 pdf(filename, width = 8, height = 10)
263 # make sure row and column names are readable
264 fontsize_row = 10 - nrow(heatmap.data) / 15
265 if (fontsize_row < 2) {
266 show_rownames = FALSE
267 }
268 else {
269 show_rownames = TRUE
270 }
271 fontsize_col = 10 - ncol(heatmap.data) / 15
272 if (fontsize_row < 2) {
273 show_colnames = FALSE
274 }
275 else {
276 show_colnames = TRUE
277 }
278 pheatmap(as.matrix(heatmap.data),
279 col = c("blue", "white", "red"),
280 legend_breaks=c(-1,0,1),
281 fontsize_row=fontsize_row,
282 fontsize_col=fontsize_col,
283 cluster_cols=TRUE,
284 main=paste("All samples, sd=", sdev, sep=""),
285 show_rownames=show_rownames,
286 show_colnames=show_colnames
287 )
288 invisible(dev.off())
289 }
290
291 loadFingerprintFromCELFile <- function(filename) {
292 ## Load the data from the provided CEL file.
293 tryCatch({ data <- ReadAffy(filenames = c(filename)) },
294 error = function(err) {
295 sink(stderr(), type = "message");
296 stop("unable to parse CEL file - ensure provided file is valid",
297 call. = FALSE)
298 })
299
300 ## Get the reported platform, if any, from the CEL file.
301 platform <- getPlatformFromArrayName(cdfName(data))
302 if (is.null(platform)) {
303 sink(stderr(), type = "message");
304 stop(sprintf("Unable to determine from platform from CEL file"));
305 };
306 if (! is.supportedPlatform(platform)) {
307 sink(stderr(), type = "message");
308 stop(sprintf("The '%s' platform is not supported by PathPrint", platform))
309 }
310
311 ## Create the fingerprint based on the CEL file data.
312 cel.fingerprint <- exprs2fingerprint(exprs(rma(data)), platform, getPlatformSpecies(platform))
313
314 ## Convert the fingerprint to a data frame.
315 cel.fingerprint <- as.data.frame(cel.fingerprint)
316
317 ## Tag the fingerprint as having been calculated on either a GSM or GSE id.
318 if (ncol(cel.fingerprint) == 1) {
319 attr(cel.fingerprint, 'fingerprintType') <- 'GSM'
320 } else {
321 attr(cel.fingerprint, 'fingerprintType') <- 'GSE'
322 }
323 attr(cel.fingerprint, 'fingerprintGEO') <- 'Unknown GEO ID'
324
325 ## Return the fingerprint
326 return(cel.fingerprint)
327 }
328
329 loadFingerprintFromExprsFile <- function(filename, platform) {
330 if (! is.supportedPlatform(platform))
331 stop(sprintf("The '%s' platform is not supported by PathPrint", platform))
332
333 ## Load the data from the provided expression set file.
334 tryCatch({ data <- read.delim(filename, sep=("\t")) },
335 error = function(err) {
336 stop(sprintf("Unable to parse expression set file: %s", err),
337 call. = FALSE)
338 })
339
340 species <- getPlatformSpecies(platform);
341 print(platform)
342 print(species)
343
344 ## Create the fingerprint based on the expression set file data.
345 tryCatch({ exprs.fingerprint <- exprs2fingerprint(as.data.frame(data), platform, getPlatformSpecies(platform)) },
346 error = function(err) {
347 print(err);
348 stop(sprintf("Expression set to fingerprint conversion failed. Please ensure the platform is correct."))
349 })
350
351 ## Convert the fingerprint to a data frame.
352 exprs.fingerprint <- as.data.frame(exprs.fingerprint)
353
354 ## Tag the fingerprint as having been calculated on either a GSM or GSE id.
355 if (ncol(exprs.fingerprint) == 1) {
356 attr(exprs.fingerprint, 'fingerprintType') <- 'GSM'
357 } else {
358 attr(exprs.fingerprint, 'fingerprintType') <- 'GSE'
359 }
360
361 ## Return the fingerprint
362 return(exprs.fingerprint)
363 }
364
365 ##
366 ## Helper functions for validating if a CEL file is supported by PathPrint
367 ## and for accessing data about the platform or species it supports.
368 ##
369
370 is.existingGeoID <- function(id) {
371 if (is.geoID(id)) {
372 if (is.geoGSM(id)) {
373 id %in% GEO.metadata.matrix[,"GSM"]
374 } else {
375 id %in% GEO.metadata.matrix[,"GSE"]
376 }
377 } else {
378 NULL
379 }
380 }
381
382 ## Check to see if the provided platform is supported by PathPrint.
383 is.supportedPlatform <- function(platform) {
384 platform %in% supportedArrays()[,1]
385 }
386
387 ## Check to see if the provided species is supported by PathPrint.
388 is.supportedSpecies <- function(species) {
389 species %in% supportedArrays()[,4]
390 }
391
392 ## Given the array name, return the corresponding species.
393 getPlatformSpecies <- function(platform) {
394 if (is.supportedPlatform(platform)) {
395 return(supportedArrays()[grep(sprintf("^%s$", platform), (supportedArrays()[,1])), 4])
396 } else {
397 return(NULL)
398 }
399 }
400
401 ## Lookup the platform using the array name.
402 getPlatformFromArrayName <- function(array) {
403 platform <- supportedArrays()[grep(sprintf("^%s$", array), (supportedArrays()[,2])), 1]
404
405 if (length(platform))
406 return(platform)
407 else
408 return(NULL)
409 }
410
411 generateSimilarExperiments <- function(data, filename) {
412 sink(filename, append=FALSE, split=FALSE)
413
414 cat("<html>\n")
415 cat("<head>\n")
416 cat("<title>Similar Experiments in GEO</title>\n")
417 cat("
418 <style type='text/css'>
419
420
421 .data-table {
422 border-collapse: collapse;
423 width: 100%;
424 }
425
426 .data-table td, th {
427 border: 1px solid lightslategrey;
428 color: #000080;
429 font-family: Verdana,Geneva,Arial,sans-serif;
430 padding: 4px;
431 font-size: 75%;
432 overflow: hidden;
433 text-overflow: ellipsis;
434 }
435
436 .data-table th {
437 font-weight: bold;
438 }
439 </style>\n");
440 cat("<head>\n")
441 cat("<body>\n")
442 cat("<table class='data-table'>\n")
443 cat("<tr>")
444 cat("<th>GSM ID</th>")
445 cat("<th>GSE ID</th>")
446 cat("<th>GPL ID</th>")
447 cat("<th>Source</th>")
448 cat("<th>distance</th>")
449 cat("<th>p-value</th>")
450 cat("</tr>\n")
451 for(i in 1:nrow(data)) {
452 cat("<tr>",
453 "<td><a href='", createGEOLink(data$GSM[i]), "' target=0>", data$GSM[i], "</a></td>",
454 "<td><a href='", createGEOLink(data$GSE[i]), "' target=0>", data$GSE[i], "</a></td>",
455 "<td><a href='", createGEOLink(data$GPL[i]), "' target=0>", data$GPL[i], "</a></td>",
456 "<td>", data$Source[i], "</td>",
457 "<td>", data$distance[i], "</td>",
458 "<td>", data$pvalue[i], "</td>",
459 "</tr>\n", sep = '')
460 }
461
462 cat("</table>\n")
463 cat("</body\n")
464 cat("</html>\n")
465
466 sink()
467 }
468
469 ##
470 ## Helper functions
471 ##
472
473 ## Create a link to the GEO ID.
474 createGEOLink <- function(id) {
475 sprintf("http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=%s", id)
476 }
477
478 ## Check if the provided ID corresponds to the GEO GSM format.
479 is.geoGSM <- function(id) {
480 length(grep("^GSM\\d+$", id, ignore.case = TRUE)) > 0
481 }
482
483 ## Check if the provided ID corresponds to the GEO GSE format.
484 is.geoGSE <- function(id) {
485 length(grep("^GSE\\d+$", id, ignore.case = TRUE)) > 0
486 }
487
488 ## Check if the provided ID corresponds to the GEO format.
489 is.geoID <- function(id) {
490 return(is.geoGSM(id) | is.geoGSE(id))
491 }
492