Mercurial > repos > iuc > limma_voom

diff limma_voom.R @ 4:a61a6e62e91f 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 6a458881c0819b75e55e64b3f494679d43bb9ee8
author iuc
date Sun, 29 Apr 2018 17:36:42 -0400
parents 38aab66ae5cb
children d8a55b5f0de0
line wrap: on
line diff
--- a/limma_voom.R	Wed Jan 31 12:45:42 2018 -0500
+++ b/limma_voom.R	Sun Apr 29 17:36:42 2018 -0400
@@ -349,7 +349,9 @@
 data <- list()
 data$counts <- counts
 if (haveAnno) {
-  data$genes <- geneanno
+  # order annotation by genes in counts (assumes gene ids are in 1st column of geneanno)
+  annoord <- geneanno[match(row.names(counts), geneanno[,1]), ]
+  data$genes <- annoord
 } else {
   data$genes <- data.frame(GeneID=row.names(counts))
 }
@@ -431,6 +433,7 @@
     # limma-trend approach
     logCPM <- cpm(data, log=TRUE, prior.count=opt$trend)
     fit <- lmFit(logCPM, design)
+    fit$genes <- data$genes
     fit <- contrasts.fit(fit, contrasts)
     if (wantRobust) {
         fit <- eBayes(fit, trend=TRUE, robust=TRUE)
@@ -459,7 +462,7 @@
 
     # Save normalised counts (log2cpm)
     if (wantNorm) {
-        write.table(logCPM, file=normOut, row.names=TRUE, sep="\t")
+        write.table(logCPM, file=normOut, row.names=TRUE, sep="\t", quote=FALSE)
         linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts.tsv"))))
     }
 } else {
@@ -510,7 +513,7 @@
      # Save normalised counts (log2cpm)
     if (wantNorm) {
         norm_counts <- data.frame(vData$genes, vData$E)
-        write.table(norm_counts, file=normOut, row.names=FALSE, sep="\t")
+        write.table(norm_counts, file=normOut, row.names=FALSE, sep="\t", quote=FALSE)
         linkData <- rbind(linkData, c((paste0(deMethod, "_", "normcounts.tsv")), (paste0(deMethod, "_", "normcounts.tsv"))))
     }
 
@@ -554,11 +557,7 @@
 
     # Write top expressions table
     top <- topTable(fit, coef=i, number=Inf, sort.by="P")
-    if (wantTrend) {
-        write.table(top, file=topOut[i], row.names=TRUE, sep="\t")
-    } else {
-        write.table(top, file=topOut[i], row.names=FALSE, sep="\t")
-    }
+    write.table(top, file=topOut[i], row.names=FALSE, sep="\t", quote=FALSE)
 
     linkName <- paste0(deMethod, "_", contrastData[i], ".tsv")
     linkAddr <- paste0(deMethod, "_", contrastData[i], ".tsv")
@@ -627,7 +626,7 @@
 
 cata("<body>\n")
 cata("<h3>Limma Analysis Output:</h3>\n")
-cata("PDF copies of JPEGS available in 'Plots' section.<br />\n")
+cata("Links to PDF copies of plots are in 'Plots' section below />\n")
 if (wantWeight) {
     HtmlImage(imageData$Link[1], imageData$Label[1], width=1000)
 } else {