Mercurial > repos > iuc > limma_voom

comparison limma_voom.R @ 7:e6a4ff41af6b 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 cf399638ebca4250bcc15f468238a9964de97b33
author iuc
date Wed, 09 May 2018 13:27:14 -0400
parents 39fa12a6d885
children 00a42b66e522
comparison
equal deleted inserted replaced
6:39fa12a6d885 7:e6a4ff41af6b
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 25 # topgenes", "G", 1, "integer" -Integer specifying no. of genes to highlight in volcano and heatmap
26 # treatOpt", "T", 0, "logical" -String specifying if TREAT function shuld be used 26 # treatOpt", "T", 0, "logical" -String specifying if TREAT function shuld be used
27 # 27 #
28 # OUT: 28 # OUT:
29 # Density Plots (if filtering)
30 # Box Plots (if normalising)
29 # MDS Plot 31 # MDS Plot
30 # Voom/SA plot 32 # Voom/SA plot
31 # MD Plot 33 # MD Plot
34 # Volcano Plot
35 # Heatmap
32 # Expression Table 36 # Expression Table
33 # HTML file linking to the ouputs 37 # HTML file linking to the ouputs
34 # Optional: 38 # Optional:
35 # Normalised counts Table 39 # Normalised counts Table
36 # RData file 40 # RData file
37 # 41 #
38 # 42 #
39 # Author: Shian Su - registertonysu@gmail.com - Jan 2014 43 # Author: Shian Su - registertonysu@gmail.com - Jan 2014
40 # Modified by: Maria Doyle - Jun 2017, Jan 2018 44 # Modified by: Maria Doyle - Jun 2017, Jan 2018, May 2018
41 45
42 # Record starting time 46 # Record starting time
43 timeStart <- as.character(Sys.time()) 47 timeStart <- as.character(Sys.time())
44 48
45 # Load all required libraries 49 # Load all required libraries
48 library(splines, quietly=TRUE, warn.conflicts=FALSE) 52 library(splines, quietly=TRUE, warn.conflicts=FALSE)
49 library(edgeR, quietly=TRUE, warn.conflicts=FALSE) 53 library(edgeR, quietly=TRUE, warn.conflicts=FALSE)
50 library(limma, quietly=TRUE, warn.conflicts=FALSE) 54 library(limma, quietly=TRUE, warn.conflicts=FALSE)
51 library(scales, quietly=TRUE, warn.conflicts=FALSE) 55 library(scales, quietly=TRUE, warn.conflicts=FALSE)
52 library(getopt, quietly=TRUE, warn.conflicts=FALSE) 56 library(getopt, quietly=TRUE, warn.conflicts=FALSE)
57 library(gplots, quietly=TRUE, warn.conflicts=FALSE)
53 58
54 ################################################################################ 59 ################################################################################
55 ### Function Declaration 60 ### Function Declaration
56 ################################################################################ 61 ################################################################################
57 # Function to sanitise contrast equations so there are no whitespaces 62 # Function to sanitise contrast equations so there are no whitespaces
167 "pAdjOpt", "d", 1, "character", 172 "pAdjOpt", "d", 1, "character",
168 "normOpt", "n", 1, "character", 173 "normOpt", "n", 1, "character",
169 "robOpt", "b", 0, "logical", 174 "robOpt", "b", 0, "logical",
170 "trend", "t", 1, "double", 175 "trend", "t", 1, "double",
171 "weightOpt", "w", 0, "logical", 176 "weightOpt", "w", 0, "logical",
172 "volhiOpt", "v", 1, "integer", 177 "topgenes", "G", 1, "integer",
173 "treatOpt", "T", 0, "logical"), 178 "treatOpt", "T", 0, "logical"),
174 byrow=TRUE, ncol=4) 179 byrow=TRUE, ncol=4)
175 opt <- getopt(spec) 180 opt <- getopt(spec)
176 181
177 182
328 mdsOutPng[i] <- makeOut(paste0("mdsplot_", names(factors)[i], ".png")) 333 mdsOutPng[i] <- makeOut(paste0("mdsplot_", names(factors)[i], ".png"))
329 } 334 }
330 335
331 mdOutPdf <- character() 336 mdOutPdf <- character()
332 volOutPdf <- character() 337 volOutPdf <- character()
338 heatOutPdf <- character()
333 mdvolOutPng <- character() 339 mdvolOutPng <- character()
334 topOut <- character() 340 topOut <- character()
335 for (i in 1:length(contrastData)) { 341 for (i in 1:length(contrastData)) {
336 mdOutPdf[i] <- makeOut(paste0("mdplot_", contrastData[i], ".pdf")) 342 mdOutPdf[i] <- makeOut(paste0("mdplot_", contrastData[i], ".pdf"))
337 volOutPdf[i] <- makeOut(paste0("volplot_", contrastData[i], ".pdf")) 343 volOutPdf[i] <- makeOut(paste0("volplot_", contrastData[i], ".pdf"))
344 heatOutPdf[i] <- makeOut(paste0("heatmap_", contrastData[i], ".pdf"))
338 mdvolOutPng[i] <- makeOut(paste0("mdvolplot_", contrastData[i], ".png")) 345 mdvolOutPng[i] <- makeOut(paste0("mdvolplot_", contrastData[i], ".png"))
339 topOut[i] <- makeOut(paste0(deMethod, "_", contrastData[i], ".tsv")) 346 topOut[i] <- makeOut(paste0(deMethod, "_", contrastData[i], ".tsv"))
340 } 347 }
341 348
342 normOut <- makeOut(paste0(deMethod, "_normcounts.tsv")) 349 normOut <- makeOut(paste0(deMethod, "_normcounts.tsv"))
689 pdf(volOutPdf[i]) 696 pdf(volOutPdf[i])
690 if (haveAnno) { 697 if (haveAnno) {
691 # labels must be in second column currently 698 # labels must be in second column currently
692 limma::volcanoplot(fit, coef=i, 699 limma::volcanoplot(fit, coef=i,
693 main=paste("Volcano Plot:", unmake.names(contrastData[i])), 700 main=paste("Volcano Plot:", unmake.names(contrastData[i])),
694 highlight=opt$volhiOpt, 701 highlight=opt$topgenes,
695 names=fit$genes[, 2]) 702 names=fit$genes[, 2])
696 } else { 703 } else {
697 limma::volcanoplot(fit, coef=i, 704 limma::volcanoplot(fit, coef=i,
698 main=paste("Volcano Plot:", unmake.names(contrastData[i])),, 705 main=paste("Volcano Plot:", unmake.names(contrastData[i])),
699 highlight=opt$volhiOpt, 706 highlight=opt$topgenes,
700 names=fit$genes$GeneID) 707 names=fit$genes$GeneID)
701 } 708 }
702 linkName <- paste0("VolcanoPlot_", contrastData[i], ".pdf") 709 linkName <- paste0("VolcanoPlot_", contrastData[i], ".pdf")
703 linkAddr <- paste0("volplot_", contrastData[i], ".pdf") 710 linkAddr <- paste0("volplot_", contrastData[i], ".pdf")
711 linkData <- rbind(linkData, c(linkName, linkAddr))
712 invisible(dev.off())
713
714 # Plot Heatmap
715 topgenes <- rownames(top[1:opt$topgenes, ])
716 if (wantTrend) {
717 topexp <- plotData[topgenes, ]
718 } else {
719 topexp <- plotData$E[topgenes, ]
720 }
721
722 pdf(heatOutPdf[i])
723 par(cex.main=0.8)
724 if (haveAnno) {
725 # labels must be in second column currently
726 heatmap.2(topexp, scale="row",
727 main=paste("Heatmap:", unmake.names(contrastData[i])),
728 col="bluered", trace="none", density.info="none",
729 margin=c(8,6), lhei=c(2,10), dendrogram="column",
730 cexRow=0.7, cexCol=0.7, srtCol=45,
731 labRow=top[topgenes, 2])
732 } else {
733 heatmap.2(topexp, scale="row",
734 main=paste("Heatmap:", unmake.names(contrastData[i])),
735 col="bluered", trace="none", density.info="none",
736 margin=c(8,6), lhei=c(2,10), dendrogram="column",
737 cexRow=0.7, cexCol=0.7, srtCol=45)
738 }
739 linkName <- paste0("Heatmap_", contrastData[i], ".pdf")
740 linkAddr <- paste0("heatmap_", contrastData[i], ".pdf")
704 linkData <- rbind(linkData, c(linkName, linkAddr)) 741 linkData <- rbind(linkData, c(linkName, linkAddr))
705 invisible(dev.off()) 742 invisible(dev.off())
706 743
707 # PNG of MD and Volcano 744 # PNG of MD and Volcano
708 png(mdvolOutPng[i], width=1200, height=600) 745 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)) 746 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", 747 limma::plotMD(fit, status=status[, i], coef=i, main="MD Plot",
711 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1), 748 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1),
712 xlab="Average Expression", ylab="logFC") 749 xlab="Average Expression", ylab="logFC")
713
714 abline(h=0, col="grey", lty=2) 750 abline(h=0, col="grey", lty=2)
715 751
716 # Volcano plots 752 # Volcano plots
717 if (haveAnno) { 753 if (haveAnno) {
718 # labels must be in second column currently 754 # labels must be in second column currently
719 limma::volcanoplot(fit, coef=i, main="Volcano Plot", 755 limma::volcanoplot(fit, coef=i, main="Volcano Plot",
720 highlight=opt$volhiOpt, 756 highlight=opt$topgenes,
721 names=fit$genes[, 2]) 757 names=fit$genes[, 2])
722 } else { 758 } else {
723 limma::volcanoplot(fit, coef=i, main="Volcano Plot", 759 limma::volcanoplot(fit, coef=i, main="Volcano Plot",
724 highlight=opt$volhiOpt, 760 highlight=opt$topgenes,
725 names=fit$genes$GeneID) 761 names=fit$genes$GeneID)
726 } 762 }
727 763
728 imgName <- paste0("MDVolPlot_", contrastData[i]) 764 imgName <- paste0("MDVolPlot_", contrastData[i])
729 imgAddr <- paste0("mdvolplot_", contrastData[i], ".png") 765 imgAddr <- paste0("mdvolplot_", contrastData[i], ".png")