changeset 30:2100480d9be5 draft

Uploaded

--- a/differential_count_models/rgedgeRpaired_nocamera.xml	Wed Aug 07 23:49:40 2013 -0400
+++ b/differential_count_models/rgedgeRpaired_nocamera.xml	Thu Aug 08 07:00:48 2013 -0400
@@ -566,8 +566,6 @@
     print("@@ using genecards substitution for urls")
     contigurls = paste0(genecards,allgenes,"\'>",allgenes,"</a>")
   }
-  print.noquote("# urls")
-  print.noquote(head(contigurls))
   print(paste("# Total low count contigs per sample = ",paste(lo,collapse=',')),quote=F) 
   cmrowsums = rowSums(workCM)
   TName=unique(group)[1]
@@ -699,12 +697,14 @@
     deSeqDatdisp = estimateDispersions(deSeqDatsizefac,fitType=DESeq_fitType)
     resDESeq = nbinomWaldTest(deSeqDatdisp, pAdjustMethod=fdrtype)
     rDESeq = as.data.frame(results(resDESeq))
+    rDESeq = cbind(Contig=rownames(workCM),rDESeq,NReads=cmrowsums)
+    srDESeq = rDESeq[order(rDESeq\$pvalue),]
+    write.table(srDESeq,file=out_DESeq2, quote=FALSE, sep="\t",row.names=F)
+    qqPlot(descr=paste(myTitle,'DESeq2 adj p qq plot'),pvector=rDESeq\$padj,outpdf='DESeq2_qqplot.pdf')
+    cat("# DESeq top 50\n")
     rDESeq = cbind(Contig=rownames(workCM),rDESeq,NReads=cmrowsums,URL=contigurls)
     srDESeq = rDESeq[order(rDESeq\$pvalue),]
-    qqPlot(descr=paste(myTitle,'DESeq2 adj p qq plot'),pvector=rDESeq\$padj,outpdf='DESeq2_qqplot.pdf')
-    cat("# DESeq top 50\n")
     print.noquote(srDESeq[1:50,])
-    write.table(srDESeq,file=out_DESeq2, quote=FALSE, sep="\t",row.names=F)
     topresults.DESeq = rDESeq[which(rDESeq\$padj < fdrthresh), ]
     DESeqcountsindex = which(allgenes %in% rownames(topresults.DESeq))
     DESeqcounts = rep(0, length(allgenes))
@@ -740,7 +740,7 @@
   }
 
   if (doVoom == T) {
-      sink('VOOM.log')
+      sink('Voom.log')
       if (doedgeR == F) {
          #### Setup DGEList object
          DGEList = DGEList(counts=workCM, group = group)
@@ -751,20 +751,22 @@
          DGEList = estimateGLMTagwiseDisp(DGEList,mydesign)
          norm.factor = DGEList\$samples\$norm.factors
          }
-      pdf("VOOM_mean_variance_plot.pdf")
+      pdf("Voom_mean_variance_plot.pdf")
       dat.voomed = voom(DGEList, mydesign, plot = TRUE, lib.size = colSums(workCM) * norm.factor)
       dev.off()
       # Use limma to fit data
       fit = lmFit(dat.voomed, mydesign)
       fit = eBayes(fit)
       rvoom = topTable(fit, coef = length(colnames(mydesign)), adj = fdrtype, n = Inf, sort="none")
-      qqPlot(descr=paste(myTitle,'VOOM-limma adj p QQ plot'),pvector=rvoom\$adj.P.Val,outpdf='VOOM_qqplot.pdf')
+      qqPlot(descr=paste(myTitle,'Voom-limma adj p QQ plot'),pvector=rvoom\$adj.P.Val,outpdf='Voom_qqplot.pdf')
       rownames(rvoom) = rownames(workCM)
+      rvoom = cbind(rvoom,NReads=cmrowsums)
+      srvoom = rvoom[order(rvoom\$P.Value),]
+      write.table(srvoom,file=out_VOOM, quote=FALSE, sep="\t",row.names=F)
       rvoom = cbind(rvoom,NReads=cmrowsums,URL=contigurls)
       srvoom = rvoom[order(rvoom\$P.Value),]
-      cat("# VOOM top 50\n")
+      cat("# Voom top 50\n")
       print(srvoom[1:50,])
-      write.table(srvoom,file=out_VOOM, quote=FALSE, sep="\t",row.names=F)
       # Use an FDR cutoff to find interesting samples for edgeR, DESeq and voom/limma
       topresults.voom = rvoom[which(rvoom\$adj.P.Val < fdrthresh), ]
       voomcountsindex = which(allgenes %in% topresults.voom\$ID)