Mercurial > repos > iuc > limma_voom

comparison test-data/out_rscript.txt @ 6:39fa12a6d885 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
planemo upload for repository https://github.com/galaxyproject/tools-iuc/tree/master/tools/limma_voom commit 60cae222b10f43f830936c19298bd723ac47e43d
author iuc
date Tue, 08 May 2018 18:12:40 -0400
parents d8a55b5f0de0
children e6a4ff41af6b
comparison
equal deleted inserted replaced
5:d8a55b5f0de0 6:39fa12a6d885
20 # pAdjOpt", "d", 1, "character" -String specifying the p-value adjustment method 20 # pAdjOpt", "d", 1, "character" -String specifying the p-value adjustment method
21 # normOpt", "n", 1, "character" -String specifying type of normalisation used 21 # normOpt", "n", 1, "character" -String specifying type of normalisation used
22 # robOpt", "b", 0, "logical" -String specifying if robust options should be used 22 # robOpt", "b", 0, "logical" -String specifying if robust options should be used
23 # trend", "t", 1, "double" -Float for prior.count if limma-trend is used instead of voom 23 # trend", "t", 1, "double" -Float for prior.count if limma-trend is used instead of voom
24 # weightOpt", "w", 0, "logical" -String specifying if voomWithQualityWeights should be used 24 # weightOpt", "w", 0, "logical" -String specifying if voomWithQualityWeights should be used
25 # volhiOpt", "v", 1, "integer" -Integer specifying number of genes to highlight in volcano plot
26 # treatOpt", "T", 0, "logical" -String specifying if TREAT function shuld be used
25 # 27 #
26 # OUT: 28 # OUT:
27 # MDS Plot 29 # MDS Plot
28 # Voom/SA plot 30 # Voom/SA plot
29 # MD Plot 31 # MD Plot
47 library(edgeR, quietly=TRUE, warn.conflicts=FALSE) 49 library(edgeR, quietly=TRUE, warn.conflicts=FALSE)
48 library(limma, quietly=TRUE, warn.conflicts=FALSE) 50 library(limma, quietly=TRUE, warn.conflicts=FALSE)
49 library(scales, quietly=TRUE, warn.conflicts=FALSE) 51 library(scales, quietly=TRUE, warn.conflicts=FALSE)
50 library(getopt, quietly=TRUE, warn.conflicts=FALSE) 52 library(getopt, quietly=TRUE, warn.conflicts=FALSE)
51 53
52 if (packageVersion("limma") < "3.20.1") {
53 stop("Please update 'limma' to version >= 3.20.1 to run this tool")
54 }
55
56 ################################################################################ 54 ################################################################################
57 ### Function Delcaration 55 ### Function Declaration
58 ################################################################################ 56 ################################################################################
59 # Function to sanitise contrast equations so there are no whitespaces 57 # Function to sanitise contrast equations so there are no whitespaces
60 # surrounding the arithmetic operators, leading or trailing whitespace 58 # surrounding the arithmetic operators, leading or trailing whitespace
61 sanitiseEquation <- function(equation) { 59 sanitiseEquation <- function(equation) {
62 equation <- gsub(" *[+] *", "+", equation) 60 equation <- gsub(" *[+] *", "+", equation)
168 "pValReq", "p", 1, "double", 166 "pValReq", "p", 1, "double",
169 "pAdjOpt", "d", 1, "character", 167 "pAdjOpt", "d", 1, "character",
170 "normOpt", "n", 1, "character", 168 "normOpt", "n", 1, "character",
171 "robOpt", "b", 0, "logical", 169 "robOpt", "b", 0, "logical",
172 "trend", "t", 1, "double", 170 "trend", "t", 1, "double",
173 "weightOpt", "w", 0, "logical"), 171 "weightOpt", "w", 0, "logical",
172 "volhiOpt", "v", 1, "integer",
173 "treatOpt", "T", 0, "logical"),
174 byrow=TRUE, ncol=4) 174 byrow=TRUE, ncol=4)
175 opt <- getopt(spec) 175 opt <- getopt(spec)
176 176
177 177
178 if (is.null(opt$matrixPath) & is.null(opt$filesPath)) { 178 if (is.null(opt$matrixPath) & is.null(opt$filesPath)) {
233 deMethod <- "limma-voom" 233 deMethod <- "limma-voom"
234 } else { 234 } else {
235 wantTrend <- TRUE 235 wantTrend <- TRUE
236 deMethod <- "limma-trend" 236 deMethod <- "limma-trend"
237 priorCount <- opt$trend 237 priorCount <- opt$trend
238 }
239
240 if (is.null(opt$treatOpt)) {
241 wantTreat <- FALSE
242 } else {
243 wantTreat <- TRUE
238 } 244 }
239 245
240 246
241 if (!is.null(opt$filesPath)) { 247 if (!is.null(opt$filesPath)) {
242 # Process the separate count files (adapted from DESeq2 wrapper) 248 # Process the separate count files (adapted from DESeq2 wrapper)
309 # Split up contrasts seperated by comma into a vector then sanitise 315 # Split up contrasts seperated by comma into a vector then sanitise
310 contrastData <- unlist(strsplit(opt$contrastData, split=",")) 316 contrastData <- unlist(strsplit(opt$contrastData, split=","))
311 contrastData <- sanitiseEquation(contrastData) 317 contrastData <- sanitiseEquation(contrastData)
312 contrastData <- gsub(" ", ".", contrastData, fixed=TRUE) 318 contrastData <- gsub(" ", ".", contrastData, fixed=TRUE)
313 319
314 320 denOutPng <- makeOut("densityplots.png")
315 mdsOutPdf <- makeOut("mdsplot_nonorm.pdf") 321 denOutPdf <- makeOut("DensityPlots.pdf")
316 mdsOutPng <- makeOut("mdsplot_nonorm.png") 322 boxOutPng <- makeOut("boxplots.png")
317 nmdsOutPdf <- makeOut("mdsplot.pdf") 323 boxOutPdf <- makeOut("BoxPlots.pdf")
318 nmdsOutPng <- makeOut("mdsplot.png") 324 mdsOutPdf <- character() # Initialise character vector
319 mdOutPdf <- character() # Initialise character vector 325 mdsOutPng <- character()
320 mdOutPng <- character() 326 for (i in 1:ncol(factors)) {
327 mdsOutPdf[i] <- makeOut(paste0("mdsplot_", names(factors)[i], ".pdf"))
328 mdsOutPng[i] <- makeOut(paste0("mdsplot_", names(factors)[i], ".png"))
329 }
330
331 mdOutPdf <- character()
332 volOutPdf <- character()
333 mdvolOutPng <- character()
321 topOut <- character() 334 topOut <- character()
322 for (i in 1:length(contrastData)) { 335 for (i in 1:length(contrastData)) {
323 mdOutPdf[i] <- makeOut(paste0("mdplot_", contrastData[i], ".pdf")) 336 mdOutPdf[i] <- makeOut(paste0("mdplot_", contrastData[i], ".pdf"))
324 mdOutPng[i] <- makeOut(paste0("mdplot_", contrastData[i], ".png")) 337 volOutPdf[i] <- makeOut(paste0("volplot_", contrastData[i], ".pdf"))
338 mdvolOutPng[i] <- makeOut(paste0("mdvolplot_", contrastData[i], ".png"))
325 topOut[i] <- makeOut(paste0(deMethod, "_", contrastData[i], ".tsv")) 339 topOut[i] <- makeOut(paste0(deMethod, "_", contrastData[i], ".tsv"))
326 } 340 }
341
327 normOut <- makeOut(paste0(deMethod, "_normcounts.tsv")) 342 normOut <- makeOut(paste0(deMethod, "_normcounts.tsv"))
328 rdaOut <- makeOut(paste0(deMethod, "_analysis.RData")) 343 rdaOut <- makeOut(paste0(deMethod, "_analysis.RData"))
329 sessionOut <- makeOut("session_info.txt") 344 sessionOut <- makeOut("session_info.txt")
330 345
331 # Initialise data for html links and images, data frame with columns Label and 346 # Initialise data for html links and images, data frame with columns Label and
380 # Creating naming data 395 # Creating naming data
381 samplenames <- colnames(data$counts) 396 samplenames <- colnames(data$counts)
382 sampleanno <- data.frame("sampleID"=samplenames, factors) 397 sampleanno <- data.frame("sampleID"=samplenames, factors)
383 398
384 399
385 # Generating the DGEList object "data" 400 # Generating the DGEList object "y"
386 print("Generating DGEList object") 401 print("Generating DGEList object")
387 data$samples <- sampleanno 402 data$samples <- sampleanno
388 data$samples$lib.size <- colSums(data$counts) 403 data$samples$lib.size <- colSums(data$counts)
389 data$samples$norm.factors <- 1 404 data$samples$norm.factors <- 1
390 row.names(data$samples) <- colnames(data$counts) 405 row.names(data$samples) <- colnames(data$counts)
391 data <- new("DGEList", data) 406 y <- new("DGEList", data)
392 407
393 print("Generating Design") 408 print("Generating Design")
394 # Name rows of factors according to their sample 409 # Name rows of factors according to their sample
395 row.names(factors) <- names(data$counts) 410 row.names(factors) <- names(data$counts)
396 factorList <- sapply(names(factors), pasteListName) 411 factorList <- sapply(names(factors), pasteListName)
404 colnames(design) <- gsub(factorList[i], "", colnames(design), fixed=TRUE) 419 colnames(design) <- gsub(factorList[i], "", colnames(design), fixed=TRUE)
405 } 420 }
406 421
407 # Calculating normalising factors 422 # Calculating normalising factors
408 print("Calculating Normalisation Factors") 423 print("Calculating Normalisation Factors")
409 data <- calcNormFactors(data, method=opt$normOpt) 424 y <- calcNormFactors(y, method=opt$normOpt)
410 425
411 # Generate contrasts information 426 # Generate contrasts information
412 print("Generating Contrasts") 427 print("Generating Contrasts")
413 contrasts <- makeContrasts(contrasts=contrastData, levels=design) 428 contrasts <- makeContrasts(contrasts=contrastData, levels=design)
414 429
415 ################################################################################ 430 ################################################################################
416 ### Data Output 431 ### Data Output
417 ################################################################################ 432 ################################################################################
433
434 # Plot Density (if filtering low counts)
435 if (filtCPM || filtSmpCount || filtTotCount) {
436 nsamples <- ncol(counts)
437
438 # PNG
439 png(denOutPng, width=1200, height=600)
440 par(mfrow=c(1,2), cex.axis=0.8)
441
442 # before filtering
443 logcpm <- cpm(counts, log=TRUE)
444 plot(density(logcpm[,1]), col=as.numeric(factors[1, 1]), lwd=2, las=2, main="", xlab="")
445 title(main="Density Plot: Raw counts", xlab="Log-cpm")
446 for (i in 2:nsamples){
447 den <- density(logcpm[,i])
448 lines(den$x, den$y, col=as.numeric(factors[i, 1]), lwd=2)
449 }
450
451 # after filtering
452 logcpm <- cpm(data$counts, log=TRUE)
453 plot(density(logcpm[,1]), col=as.numeric(factors[1, 1]), lwd=2, las=2, main="", xlab="")
454 title(main="Density Plot: Filtered counts", xlab="Log-cpm")
455 for (i in 2:nsamples){
456 den <- density(logcpm[,i])
457 lines(den$x, den$y, col=as.numeric(factors[i, 1]), lwd=2)
458 }
459 legend("topright", samplenames, text.col=as.numeric(factors[, 1]), bty="n")
460 imgName <- "Densityplots.png"
461 imgAddr <- "densityplots.png"
462 imageData <- rbind(imageData, data.frame(Label=imgName, Link=imgAddr, stringsAsFactors=FALSE))
463 invisible(dev.off())
464
465 # PDF
466 pdf(denOutPdf, width=14)
467 par(mfrow=c(1,2), cex.axis=0.8)
468 logcpm <- cpm(counts, log=TRUE)
469 plot(density(logcpm[,1]), col=as.numeric(factors[1, 1]), lwd=2, las=2, main="", xlab="")
470 title(main="Density Plot: Raw counts", xlab="Log-cpm")
471 for (i in 2:nsamples){
472 den <- density(logcpm[,i])
473 lines(den$x, den$y, col=as.numeric(factors[i, 1]), lwd=2)
474 }
475 logcpm <- cpm(data$counts, log=TRUE)
476 plot(density(logcpm[,1]), col=as.numeric(factors[1, 1]), lwd=2, las=2, main="", xlab="")
477 title(main="Density Plot: Filtered counts", xlab="Log-cpm")
478 for (i in 2:nsamples){
479 den <- density(logcpm[,i])
480 lines(den$x, den$y, col=as.numeric(factors[i, 1]), lwd=2)
481 }
482 legend("topright", samplenames, text.col=as.numeric(factors[, 1]), bty="n")
483 linkName <- "DensityPlots.pdf"
484 linkAddr <- "densityplots.pdf"
485 linkData <- rbind(linkData, data.frame(Label=linkName, Link=linkAddr, stringsAsFactors=FALSE))
486 invisible(dev.off())
487 }
488
489 # Plot Box plots (before and after normalisation)
490 if (opt$normOpt != "none") {
491 png(boxOutPng, width=1200, height=600)
492 par(mfrow=c(1,2), cex.axis=0.8)
493 logcpm <- cpm(y$counts, log=TRUE)
494 boxplot(logcpm, las=2, col=as.numeric(factors[, 1]))
495 abline(h=median(logcpm), col=4)
496 title(main="Box Plot: Unnormalised counts", ylab="Log-cpm")
497 lcpm <- cpm(y, log=TRUE)
498 boxplot(lcpm, las=2, col=as.numeric(factors[, 1]))
499 abline(h=median(lcpm), col=4)
500 title(main="Box Plot: Normalised counts", ylab="Log-cpm")
501 imgName <- "Boxplots.png"
502 imgAddr <- "boxplots.png"
503 imageData <- rbind(imageData, data.frame(Label=imgName, Link=imgAddr, stringsAsFactors=FALSE))
504 invisible(dev.off())
505
506 pdf(boxOutPdf, width=14)
507 par(mfrow=c(1,2), cex.axis=0.8)
508 logcpm <- cpm(y$counts, log=TRUE)
509 boxplot(logcpm, las=2, col=as.numeric(factors[, 1]))
510 abline(h=median(logcpm), col=4)
511 title(main="Box Plot: Unnormalised counts", ylab="Log-cpm")
512 lcpm <- cpm(y, log=TRUE)
513 boxplot(lcpm, las=2, col=as.numeric(factors[, 1]))
514 abline(h=median(lcpm), col=4)
515 title(main="Box Plot: Normalised counts", ylab="Log-cpm")
516 linkName <- "BoxPlots.pdf"
517 linkAddr <- "boxplots.pdf"
518 linkData <- rbind(linkData, data.frame(Label=linkName, Link=linkAddr, stringsAsFactors=FALSE))
519 invisible(dev.off())
520 }
521
418 # Plot MDS 522 # Plot MDS
419 print("Generating MDS plot") 523 print("Generating MDS plot")
420 labels <- names(counts) 524 labels <- names(counts)
421 png(mdsOutPng, width=600, height=600) 525
422 # Currently only using a single factor 526 for (i in 1:ncol(factors)) {
423 plotMDS(data, labels=labels, col=as.numeric(factors[, 1]), cex=0.8, main="MDS Plot (unnormalised)") 527 png(mdsOutPng[i], width=600, height=600)
424 imageData[1, ] <- c("MDS Plot (unnormalised)", "mdsplot_nonorm.png") 528 plotMDS(y, labels=labels, col=as.numeric(factors[, i]), cex=0.8, main=paste("MDS Plot:", names(factors)[i]))
425 invisible(dev.off()) 529 imgName <- paste0("MDSPlot_", names(factors)[i], ".png")
426 530 imgAddr <- paste0("mdsplot_", names(factors)[i], ".png")
427 pdf(mdsOutPdf) 531 imageData <- rbind(imageData, data.frame(Label=imgName, Link=imgAddr, stringsAsFactors=FALSE))
428 plotMDS(data, labels=labels, cex=0.5) 532 invisible(dev.off())
429 linkData[1, ] <- c("MDS Plot (unnormalised).pdf", "mdsplot_nonorm.pdf") 533
430 invisible(dev.off()) 534 pdf(mdsOutPdf[i])
535 plotMDS(y, labels=labels, col=as.numeric(factors[, i]), cex=0.8, main=paste("MDS Plot:", names(factors)[i]))
536 linkName <- paste0("MDSPlot_", names(factors)[i], ".pdf")
537 linkAddr <- paste0("mdsplot_", names(factors)[i], ".pdf")
538 linkData <- rbind(linkData, data.frame(Label=linkName, Link=linkAddr, stringsAsFactors=FALSE))
539 invisible(dev.off())
540 }
431 541
432 if (wantTrend) { 542 if (wantTrend) {
433 # limma-trend approach 543 # limma-trend approach
434 logCPM <- cpm(data, log=TRUE, prior.count=opt$trend) 544 logCPM <- cpm(y, log=TRUE, prior.count=opt$trend)
435 fit <- lmFit(logCPM, design) 545 fit <- lmFit(logCPM, design)
436 fit$genes <- data$genes 546 fit$genes <- y$genes
437 fit <- contrasts.fit(fit, contrasts) 547 fit <- contrasts.fit(fit, contrasts)
438 if (wantRobust) { 548 if (wantRobust) {
439 fit <- eBayes(fit, trend=TRUE, robust=TRUE) 549 fit <- eBayes(fit, trend=TRUE, robust=TRUE)
440 } else { 550 } else {
441 fit <- eBayes(fit, trend=TRUE, robust=FALSE) 551 fit <- eBayes(fit, trend=TRUE, robust=FALSE)
444 saOutPng <- makeOut("saplot.png") 554 saOutPng <- makeOut("saplot.png")
445 saOutPdf <- makeOut("saplot.pdf") 555 saOutPdf <- makeOut("saplot.pdf")
446 556
447 png(saOutPng, width=600, height=600) 557 png(saOutPng, width=600, height=600)
448 plotSA(fit, main="SA Plot") 558 plotSA(fit, main="SA Plot")
449 imgName <- "SA Plot.png" 559 imgName <- "SAPlot.png"
450 imgAddr <- "saplot.png" 560 imgAddr <- "saplot.png"
451 imageData <- rbind(imageData, c(imgName, imgAddr)) 561 imageData <- rbind(imageData, c(imgName, imgAddr))
452 invisible(dev.off()) 562 invisible(dev.off())
453 563
454 pdf(saOutPdf, width=14) 564 pdf(saOutPdf, width=14)
455 plotSA(fit, main="SA Plot") 565 plotSA(fit, main="SA Plot")
456 linkName <- paste0("SA Plot.pdf") 566 linkName <- "SAPlot.pdf"
457 linkAddr <- paste0("saplot.pdf") 567 linkAddr <- "saplot.pdf"
458 linkData <- rbind(linkData, c(linkName, linkAddr)) 568 linkData <- rbind(linkData, c(linkName, linkAddr))
459 invisible(dev.off()) 569 invisible(dev.off())
460 570
461 plotData <- logCPM 571 plotData <- logCPM
462 572
471 voomOutPng <- makeOut("voomplot.png") 581 voomOutPng <- makeOut("voomplot.png")
472 582
473 if (wantWeight) { 583 if (wantWeight) {
474 # Creating voom data object and plot 584 # Creating voom data object and plot
475 png(voomOutPng, width=1000, height=600) 585 png(voomOutPng, width=1000, height=600)
476 vData <- voomWithQualityWeights(data, design=design, plot=TRUE) 586 vData <- voomWithQualityWeights(y, design=design, plot=TRUE)
477 imgName <- "Voom Plot.png" 587 imgName <- "VoomPlot.png"
478 imgAddr <- "voomplot.png" 588 imgAddr <- "voomplot.png"
479 imageData <- rbind(imageData, c(imgName, imgAddr)) 589 imageData <- rbind(imageData, c(imgName, imgAddr))
480 invisible(dev.off()) 590 invisible(dev.off())
481 591
482 pdf(voomOutPdf, width=14) 592 pdf(voomOutPdf, width=14)
483 vData <- voomWithQualityWeights(data, design=design, plot=TRUE) 593 vData <- voomWithQualityWeights(y, design=design, plot=TRUE)
484 linkName <- paste0("Voom Plot.pdf") 594 linkName <- "VoomPlot.pdf"
485 linkAddr <- paste0("voomplot.pdf") 595 linkAddr <- "voomplot.pdf"
486 linkData <- rbind(linkData, c(linkName, linkAddr)) 596 linkData <- rbind(linkData, c(linkName, linkAddr))
487 invisible(dev.off()) 597 invisible(dev.off())
488 598
489 # Generating fit data and top table with weights 599 # Generating fit data and top table with weights
490 wts <- vData$weights 600 wts <- vData$weights
491 voomFit <- lmFit(vData, design, weights=wts) 601 voomFit <- lmFit(vData, design, weights=wts)
492 602
493 } else { 603 } else {
494 # Creating voom data object and plot 604 # Creating voom data object and plot
495 png(voomOutPng, width=600, height=600) 605 png(voomOutPng, width=600, height=600)
496 vData <- voom(data, design=design, plot=TRUE) 606 vData <- voom(y, design=design, plot=TRUE)
497 imgName <- "Voom Plot" 607 imgName <- "VoomPlot"
498 imgAddr <- "voomplot.png" 608 imgAddr <- "voomplot.png"
499 imageData <- rbind(imageData, c(imgName, imgAddr)) 609 imageData <- rbind(imageData, c(imgName, imgAddr))
500 invisible(dev.off()) 610 invisible(dev.off())
501 611
502 pdf(voomOutPdf) 612 pdf(voomOutPdf)
503 vData <- voom(data, design=design, plot=TRUE) 613 vData <- voom(y, design=design, plot=TRUE)
504 linkName <- paste0("Voom Plot.pdf") 614 linkName <- "VoomPlot.pdf"
505 linkAddr <- paste0("voomplot.pdf") 615 linkAddr <- "voomplot.pdf"
506 linkData <- rbind(linkData, c(linkName, linkAddr)) 616 linkData <- rbind(linkData, c(linkName, linkAddr))
507 invisible(dev.off()) 617 invisible(dev.off())
508 618
509 # Generate voom fit 619 # Generate voom fit
510 voomFit <- lmFit(vData, design) 620 voomFit <- lmFit(vData, design)
525 fit <- eBayes(voomFit, robust=FALSE) 635 fit <- eBayes(voomFit, robust=FALSE)
526 } 636 }
527 plotData <- vData 637 plotData <- vData
528 } 638 }
529 639
530 print("Generating normalised MDS plot")
531 png(nmdsOutPng, width=600, height=600)
532 # Currently only using a single factor
533 plotMDS(plotData, labels=labels, col=as.numeric(factors[, 1]), cex=0.8, main="MDS Plot (normalised)")
534 imgName <- "MDS Plot (normalised)"
535 imgAddr <- "mdsplot.png"
536 imageData <- rbind(imageData, c(imgName, imgAddr))
537 invisible(dev.off())
538
539 pdf(nmdsOutPdf)
540 plotMDS(plotData, labels=labels, cex=0.5)
541 linkName <- paste0("MDS Plot (normalised).pdf")
542 linkAddr <- paste0("mdsplot.pdf")
543 linkData <- rbind(linkData, c(linkName, linkAddr))
544 invisible(dev.off())
545
546 640
547 print("Generating DE results") 641 print("Generating DE results")
642
643 if (wantTreat) {
644 print("Applying TREAT method")
645 if (wantRobust) {
646 fit <- treat(fit, lfc=opt$lfcReq, robust=TRUE)
647 } else {
648 fit <- treat(fit, lfc=opt$lfcReq, robust=FALSE)
649 }
650 }
651
548 status = decideTests(fit, adjust.method=opt$pAdjOpt, p.value=opt$pValReq, 652 status = decideTests(fit, adjust.method=opt$pAdjOpt, p.value=opt$pValReq,
549 lfc=opt$lfcReq) 653 lfc=opt$lfcReq)
550 sumStatus <- summary(status) 654 sumStatus <- summary(status)
551 655
552 for (i in 1:length(contrastData)) { 656 for (i in 1:length(contrastData)) {
554 upCount[i] <- sumStatus["Up", i] 658 upCount[i] <- sumStatus["Up", i]
555 downCount[i] <- sumStatus["Down", i] 659 downCount[i] <- sumStatus["Down", i]
556 flatCount[i] <- sumStatus["NotSig", i] 660 flatCount[i] <- sumStatus["NotSig", i]
557 661
558 # Write top expressions table 662 # Write top expressions table
559 top <- topTable(fit, coef=i, number=Inf, sort.by="P") 663 if (wantTreat) {
664 top <- topTreat(fit, coef=i, number=Inf, sort.by="P")
665 } else{
666 top <- topTable(fit, coef=i, number=Inf, sort.by="P")
667 }
560 write.table(top, file=topOut[i], row.names=FALSE, sep="\t", quote=FALSE) 668 write.table(top, file=topOut[i], row.names=FALSE, sep="\t", quote=FALSE)
561 669
562 linkName <- paste0(deMethod, "_", contrastData[i], ".tsv") 670 linkName <- paste0(deMethod, "_", contrastData[i], ".tsv")
563 linkAddr <- paste0(deMethod, "_", contrastData[i], ".tsv") 671 linkAddr <- paste0(deMethod, "_", contrastData[i], ".tsv")
564 linkData <- rbind(linkData, c(linkName, linkAddr)) 672 linkData <- rbind(linkData, c(linkName, linkAddr))
570 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1), 678 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1),
571 xlab="Average Expression", ylab="logFC") 679 xlab="Average Expression", ylab="logFC")
572 680
573 abline(h=0, col="grey", lty=2) 681 abline(h=0, col="grey", lty=2)
574 682
575 linkName <- paste0("MD Plot_", contrastData[i], " (.pdf)") 683 linkName <- paste0("MDPlot_", contrastData[i], ".pdf")
576 linkAddr <- paste0("mdplot_", contrastData[i], ".pdf") 684 linkAddr <- paste0("mdplot_", contrastData[i], ".pdf")
577 linkData <- rbind(linkData, c(linkName, linkAddr)) 685 linkData <- rbind(linkData, c(linkName, linkAddr))
578 invisible(dev.off()) 686 invisible(dev.off())
579 687
580 png(mdOutPng[i], height=600, width=600) 688 # Plot Volcano
581 limma::plotMD(fit, status=status[, i], coef=i, 689 pdf(volOutPdf[i])
582 main=paste("MD Plot:", unmake.names(contrastData[i])), 690 if (haveAnno) {
583 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1), 691 # labels must be in second column currently
584 xlab="Average Expression", ylab="logFC") 692 limma::volcanoplot(fit, coef=i,
693 main=paste("Volcano Plot:", unmake.names(contrastData[i])),
694 highlight=opt$volhiOpt,
695 names=fit$genes[, 2])
696 } else {
697 limma::volcanoplot(fit, coef=i,
698 main=paste("Volcano Plot:", unmake.names(contrastData[i])),,
699 highlight=opt$volhiOpt,
700 names=fit$genes$GeneID)
701 }
702 linkName <- paste0("VolcanoPlot_", contrastData[i], ".pdf")
703 linkAddr <- paste0("volplot_", contrastData[i], ".pdf")
704 linkData <- rbind(linkData, c(linkName, linkAddr))
705 invisible(dev.off())
706
707 # PNG of MD and Volcano
708 png(mdvolOutPng[i], width=1200, height=600)
709 par(mfrow=c(1, 2), mar=c(5,4,2,2)+0.1, oma=c(0,0,3,0))
710 limma::plotMD(fit, status=status[, i], coef=i, main="MD Plot",
711 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1),
712 xlab="Average Expression", ylab="logFC")
585 713
586 abline(h=0, col="grey", lty=2) 714 abline(h=0, col="grey", lty=2)
587 715
588 imgName <- paste0("MD Plot_", contrastData[i]) 716 # Volcano plots
589 imgAddr <- paste0("mdplot_", contrastData[i], ".png") 717 if (haveAnno) {
718 # labels must be in second column currently
719 limma::volcanoplot(fit, coef=i, main="Volcano Plot",
720 highlight=opt$volhiOpt,
721 names=fit$genes[, 2])
722 } else {
723 limma::volcanoplot(fit, coef=i, main="Volcano Plot",
724 highlight=opt$volhiOpt,
725 names=fit$genes$GeneID)
726 }
727
728 imgName <- paste0("MDVolPlot_", contrastData[i])
729 imgAddr <- paste0("mdvolplot_", contrastData[i], ".png")
590 imageData <- rbind(imageData, c(imgName, imgAddr)) 730 imageData <- rbind(imageData, c(imgName, imgAddr))
731 title(paste0("Contrast: ", unmake.names(contrastData[i])), outer=TRUE, cex.main=1.5)
591 invisible(dev.off()) 732 invisible(dev.off())
592 } 733 }
593 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount) 734 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount)
594 row.names(sigDiff) <- contrastData 735 row.names(sigDiff) <- contrastData
595 736
596 # Save relevant items as rda object 737 # Save relevant items as rda object
597 if (wantRda) { 738 if (wantRda) {
598 print("Saving RData") 739 print("Saving RData")
599 if (wantWeight) { 740 if (wantWeight) {
600 save(data, status, plotData, labels, factors, wts, fit, top, contrasts, 741 save(counts, data, y, status, plotData, labels, factors, wts, fit, top, contrasts, design,
601 design,
602 file=rdaOut, ascii=TRUE) 742 file=rdaOut, ascii=TRUE)
603 } else { 743 } else {
604 save(data, status, plotData, labels, factors, fit, top, contrasts, design, 744 save(counts, data, y, status, plotData, labels, factors, fit, top, contrasts, design,
605 file=rdaOut, ascii=TRUE) 745 file=rdaOut, ascii=TRUE)
606 } 746 }
607 linkData <- rbind(linkData, c((paste0(deMethod, "_analysis.RData")), (paste0(deMethod, "_analysis.RData")))) 747 linkData <- rbind(linkData, c((paste0(deMethod, "_analysis.RData")), (paste0(deMethod, "_analysis.RData"))))
608 } 748 }
609 749
624 764
625 cata("<html>\n") 765 cata("<html>\n")
626 766
627 cata("<body>\n") 767 cata("<body>\n")
628 cata("<h3>Limma Analysis Output:</h3>\n") 768 cata("<h3>Limma Analysis Output:</h3>\n")
629 cata("Links to PDF copies of plots are in 'Plots' section below />\n") 769 cata("Links to PDF copies of plots are in 'Plots' section below <br />\n")
630 if (wantWeight) { 770
631 HtmlImage(imageData$Link[1], imageData$Label[1], width=1000) 771 for (i in 1:nrow(imageData)) {
632 } else { 772 if (grepl("density|box|mdvol", imageData$Link[i])) {
633 HtmlImage(imageData$Link[1], imageData$Label[1]) 773 HtmlImage(imageData$Link[i], imageData$Label[i], width=1200)
634 } 774 } else if (wantWeight) {
635 775 HtmlImage(imageData$Link[i], imageData$Label[i], width=1000)
636 for (i in 2:nrow(imageData)) { 776 } else {
637 HtmlImage(imageData$Link[i], imageData$Label[i]) 777 HtmlImage(imageData$Link[i], imageData$Label[i])
778 }
638 } 779 }
639 780
640 cata("<h4>Differential Expression Counts:</h4>\n") 781 cata("<h4>Differential Expression Counts:</h4>\n")
641 782
642 cata("<table border=\"1\" cellpadding=\"4\">\n") 783 cata("<table border=\"1\" cellpadding=\"4\">\n")
717 if (wantWeight) { 858 if (wantWeight) {
718 ListItem("Weights were applied to samples.") 859 ListItem("Weights were applied to samples.")
719 } else { 860 } else {
720 ListItem("Weights were not applied to samples.") 861 ListItem("Weights were not applied to samples.")
721 } 862 }
863 if (wantTreat) {
864 ListItem(paste("Testing significance relative to a fold-change threshold (TREAT) was performed using a threshold of log2 =", opt$lfcReq, "at FDR of", opt$pValReq, "."))
865 }
722 if (wantRobust) { 866 if (wantRobust) {
723 ListItem("eBayes was used with robust settings (robust=TRUE).") 867 if (wantTreat) {
868 ListItem("TREAT was used with robust settings (robust=TRUE).")
869 } else {
870 ListItem("eBayes was used with robust settings (robust=TRUE).")
871 }
724 } 872 }
725 if (opt$pAdjOpt!="none") { 873 if (opt$pAdjOpt!="none") {
726 if (opt$pAdjOpt=="BH" || opt$pAdjOpt=="BY") { 874 if (opt$pAdjOpt=="BH" || opt$pAdjOpt=="BY") {
727 tempStr <- paste0("MD Plot highlighted genes are significant at FDR ", 875 tempStr <- paste0("MD Plot highlighted genes are significant at FDR ",
728 "of ", opt$pValReq," and exhibit log2-fold-change of at ", 876 "of ", opt$pValReq," and exhibit log2-fold-change of at ",