Mercurial > repos > iuc > limma_voom

comparison limma_voom.R @ 5:d8a55b5f0de0 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 5f0db9544002b02d3ccffb8c983b5c8baef6bacd
author iuc
date Sat, 05 May 2018 17:55:13 -0400
parents a61a6e62e91f
children 39fa12a6d885
comparison
equal deleted inserted replaced
4:a61a6e62e91f 5:d8a55b5f0de0
314 314
315 mdsOutPdf <- makeOut("mdsplot_nonorm.pdf") 315 mdsOutPdf <- makeOut("mdsplot_nonorm.pdf")
316 mdsOutPng <- makeOut("mdsplot_nonorm.png") 316 mdsOutPng <- makeOut("mdsplot_nonorm.png")
317 nmdsOutPdf <- makeOut("mdsplot.pdf") 317 nmdsOutPdf <- makeOut("mdsplot.pdf")
318 nmdsOutPng <- makeOut("mdsplot.png") 318 nmdsOutPng <- makeOut("mdsplot.png")
319 maOutPdf <- character() # Initialise character vector 319 mdOutPdf <- character() # Initialise character vector
320 maOutPng <- character() 320 mdOutPng <- character()
321 topOut <- character() 321 topOut <- character()
322 for (i in 1:length(contrastData)) { 322 for (i in 1:length(contrastData)) {
323 maOutPdf[i] <- makeOut(paste0("maplot_", contrastData[i], ".pdf")) 323 mdOutPdf[i] <- makeOut(paste0("mdplot_", contrastData[i], ".pdf"))
324 maOutPng[i] <- makeOut(paste0("maplot_", contrastData[i], ".png")) 324 mdOutPng[i] <- makeOut(paste0("mdplot_", contrastData[i], ".png"))
325 topOut[i] <- makeOut(paste0(deMethod, "_", contrastData[i], ".tsv")) 325 topOut[i] <- makeOut(paste0(deMethod, "_", contrastData[i], ".tsv"))
326 } 326 }
327 normOut <- makeOut(paste0(deMethod, "_normcounts.tsv")) 327 normOut <- makeOut(paste0(deMethod, "_normcounts.tsv"))
328 rdaOut <- makeOut(paste0(deMethod, "_analysis.RData")) 328 rdaOut <- makeOut(paste0(deMethod, "_analysis.RData"))
329 sessionOut <- makeOut("session_info.txt") 329 sessionOut <- makeOut("session_info.txt")
561 561
562 linkName <- paste0(deMethod, "_", contrastData[i], ".tsv") 562 linkName <- paste0(deMethod, "_", contrastData[i], ".tsv")
563 linkAddr <- paste0(deMethod, "_", contrastData[i], ".tsv") 563 linkAddr <- paste0(deMethod, "_", contrastData[i], ".tsv")
564 linkData <- rbind(linkData, c(linkName, linkAddr)) 564 linkData <- rbind(linkData, c(linkName, linkAddr))
565 565
566 # Plot MA (log ratios vs mean average) using limma package on weighted 566 # Plot MD (log ratios vs mean average) using limma package on weighted
567 pdf(maOutPdf[i]) 567 pdf(mdOutPdf[i])
568 limma::plotMD(fit, status=status, coef=i, 568 limma::plotMD(fit, status=status[, i], coef=i,
569 main=paste("MA Plot:", unmake.names(contrastData[i])), 569 main=paste("MD Plot:", unmake.names(contrastData[i])),
570 col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"), 570 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1),
571 xlab="Average Expression", ylab="logFC") 571 xlab="Average Expression", ylab="logFC")
572 572
573 abline(h=0, col="grey", lty=2) 573 abline(h=0, col="grey", lty=2)
574 574
575 linkName <- paste0("MA Plot_", contrastData[i], " (.pdf)") 575 linkName <- paste0("MD Plot_", contrastData[i], " (.pdf)")
576 linkAddr <- paste0("maplot_", contrastData[i], ".pdf") 576 linkAddr <- paste0("mdplot_", contrastData[i], ".pdf")
577 linkData <- rbind(linkData, c(linkName, linkAddr)) 577 linkData <- rbind(linkData, c(linkName, linkAddr))
578 invisible(dev.off()) 578 invisible(dev.off())
579 579
580 png(maOutPng[i], height=600, width=600) 580 png(mdOutPng[i], height=600, width=600)
581 limma::plotMD(fit, status=status, coef=i, 581 limma::plotMD(fit, status=status[, i], coef=i,
582 main=paste("MA Plot:", unmake.names(contrastData[i])), 582 main=paste("MD Plot:", unmake.names(contrastData[i])),
583 col=alpha(c("firebrick", "blue"), 0.4), values=c("1", "-1"), 583 hl.col=alpha(c("firebrick", "blue"), 0.4), values=c(1, -1),
584 xlab="Average Expression", ylab="logFC") 584 xlab="Average Expression", ylab="logFC")
585 585
586 abline(h=0, col="grey", lty=2) 586 abline(h=0, col="grey", lty=2)
587 587
588 imgName <- paste0("MA Plot_", contrastData[i]) 588 imgName <- paste0("MD Plot_", contrastData[i])
589 imgAddr <- paste0("maplot_", contrastData[i], ".png") 589 imgAddr <- paste0("mdplot_", contrastData[i], ".png")
590 imageData <- rbind(imageData, c(imgName, imgAddr)) 590 imageData <- rbind(imageData, c(imgName, imgAddr))
591 invisible(dev.off()) 591 invisible(dev.off())
592 } 592 }
593 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount) 593 sigDiff <- data.frame(Up=upCount, Flat=flatCount, Down=downCount)
594 row.names(sigDiff) <- contrastData 594 row.names(sigDiff) <- contrastData
687 cata("<h4>Additional Information</h4>\n") 687 cata("<h4>Additional Information</h4>\n")
688 cata("<ul>\n") 688 cata("<ul>\n")
689 689
690 if (filtCPM || filtSmpCount || filtTotCount) { 690 if (filtCPM || filtSmpCount || filtTotCount) {
691 if (filtCPM) { 691 if (filtCPM) {
692 tempStr <- paste("Genes without more than", opt$cmpReq, 692 tempStr <- paste("Genes without more than", opt$cpmReq,
693 "CPM in at least", opt$sampleReq, "samples are insignificant", 693 "CPM in at least", opt$sampleReq, "samples are insignificant",
694 "and filtered out.") 694 "and filtered out.")
695 } else if (filtSmpCount) { 695 } else if (filtSmpCount) {
696 tempStr <- paste("Genes without more than", opt$cntReq, 696 tempStr <- paste("Genes without more than", opt$cntReq,
697 "counts in at least", opt$sampleReq, "samples are insignificant", 697 "counts in at least", opt$sampleReq, "samples are insignificant",
722 if (wantRobust) { 722 if (wantRobust) {
723 ListItem("eBayes was used with robust settings (robust=TRUE).") 723 ListItem("eBayes was used with robust settings (robust=TRUE).")
724 } 724 }
725 if (opt$pAdjOpt!="none") { 725 if (opt$pAdjOpt!="none") {
726 if (opt$pAdjOpt=="BH" || opt$pAdjOpt=="BY") { 726 if (opt$pAdjOpt=="BH" || opt$pAdjOpt=="BY") {
727 tempStr <- paste0("MA-Plot highlighted genes are significant at FDR ", 727 tempStr <- paste0("MD Plot highlighted genes are significant at FDR ",
728 "of ", opt$pValReq," and exhibit log2-fold-change of at ", 728 "of ", opt$pValReq," and exhibit log2-fold-change of at ",
729 "least ", opt$lfcReq, ".") 729 "least ", opt$lfcReq, ".")
730 ListItem(tempStr) 730 ListItem(tempStr)
731 } else if (opt$pAdjOpt=="holm") { 731 } else if (opt$pAdjOpt=="holm") {
732 tempStr <- paste0("MA-Plot highlighted genes are significant at adjusted ", 732 tempStr <- paste0("MD Plot highlighted genes are significant at adjusted ",
733 "p-value of ", opt$pValReq," by the Holm(1979) ", 733 "p-value of ", opt$pValReq," by the Holm(1979) ",
734 "method, and exhibit log2-fold-change of at least ", 734 "method, and exhibit log2-fold-change of at least ",
735 opt$lfcReq, ".") 735 opt$lfcReq, ".")
736 ListItem(tempStr) 736 ListItem(tempStr)
737 } 737 }
738 } else { 738 } else {
739 tempStr <- paste0("MA-Plot highlighted genes are significant at p-value ", 739 tempStr <- paste0("MD Plot highlighted genes are significant at p-value ",
740 "of ", opt$pValReq," and exhibit log2-fold-change of at ", 740 "of ", opt$pValReq," and exhibit log2-fold-change of at ",
741 "least ", opt$lfcReq, ".") 741 "least ", opt$lfcReq, ".")
742 ListItem(tempStr) 742 ListItem(tempStr)
743 } 743 }
744 cata("</ul>\n") 744 cata("</ul>\n")