Mercurial > repos > yhoogstrate > edger_with_design_matrix

changeset 15:db9eb4b6f778 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author yhoogstrate
date Thu, 09 Jan 2014 08:51:59 -0500
parents 15fd1a6798e0
children fba5577122a8
files edgeR_DGE.xml
diffstat 1 files changed, 20 insertions(+), 21 deletions(-) [+]
line wrap: on
line diff
--- a/edgeR_DGE.xml	Thu Jan 09 08:29:43 2014 -0500
+++ b/edgeR_DGE.xml	Thu Jan 09 08:51:59 2014 -0500
@@ -8,7 +8,7 @@
 			http://www.cheetahtemplate.org/docs/users_guide_html_multipage/contents.html
 		-->
 		
-		R CMD BATCH --vanilla --slave '--args
+		R --vanilla --slave -f $R_script '--args
 			$design_matrix
 			$contrast
 			
@@ -22,7 +22,7 @@
 			$output_BCVplot
 			$output_MAplot
 			smearPlot '
-			$R_script $output_R
+			> $output_R
 	</command>
 	
 	<inputs>
@@ -52,7 +52,7 @@
 
 output_1            = args[3]
 output_2            = args[4]
-output_3            = args[5]		##FPKM file - to be implemented
+output_3            = args[5]		##FPKM file - yet to be implemented
 output_4            = args[6]
 
 QC                  = nchar(args[7]) > 0
@@ -75,15 +75,14 @@
 read_counts = read.delim(as.character(raw_data[1,1]),header=has_header,stringsAsFactors=F,row.names=1)[1]
 
 for(i in 2:length(raw_data[,1])) {
-  print("parsing counts from:")
-  print(raw_data[i,1])
+  write("parsing counts from:",stdout())
+  write(raw_data[i,1],stdout())
   
   header = read.delim(as.character(raw_data[i,1]),header=F,stringsAsFactors=F,row.names=1,nrows=1)
   has_header = (class(header[1,1]) == "character")
   table = read.delim(as.character(raw_data[i,1]),header=has_header,stringsAsFactors=F,row.names=1)[1]
   
   read_counts = cbind(read_counts,table)
-  print(i)
 }
 
 colnames(read_counts) = as.character(raw_data[,2])
@@ -122,20 +121,20 @@
 
 
 
-  print("Calculating normalization factors...")
+  write("Calculating normalization factors...",stdout())
   dge = calcNormFactors(dge)
-  print("Estimating common dispersion...")
+  write("Estimating common dispersion...",stdout())
   dge = estimateGLMCommonDisp(dge,design)
-  print("Estimating trended dispersion...")
+  write("Estimating trended dispersion...",stdout())
   dge = estimateGLMTrendedDisp(dge,design)
-  print("Estimating tagwise dispersion...")
+  write("Estimating tagwise dispersion...",stdout())
   dge = estimateGLMTagwiseDisp(dge,design)
 
 
 
 
-  if (QC == TRUE) {
-    print("Creating QC plots...")
+  if(QC == TRUE) {
+    write("Creating QC plots...",stdout())
     #### MDS Plot
     pdf(output_5)
     plotMDS(dge, main="edgeR MDS Plot")
@@ -148,24 +147,24 @@
 
 
 
-  print("Fitting GLM...")
+  write("Fitting GLM...",stdout())
   fit   = glmFit(dge,design)
 
-  print(paste("Performing likelihood ratio test: ",contrast,sep=""))
+  write(paste("Performing likelihood ratio test: ",contrast,sep=""),stdout())
   cont &lt;- c(contrast)
   cont &lt;- makeContrasts(contrasts=cont, levels=design)
 
   lrt &lt;- glmLRT(fit, contrast=cont[,1])
-  print(paste("Exporting to file: ",output_1,sep=""))
+  write(paste("Exporting to file: ",output_1,sep=""),stdout())
   write.table(file=output_1,topTags(lrt,n=nrow(read_counts))\$table,sep="\t",row.names=T)
   write.table(file=output_2,cpm(dge,normalized.lib.sizes=TRUE),sep="\t")
   ## todo EXPORT FPKM
   write.table(file=output_4,dge\$counts,sep="\t")
-
-
-
-  if (QC == TRUE) {
-    print("Creating MA plots...")
+  
+  
+  
+  if(QC == TRUE) {
+    write("Creating MA plots...",stdout())
     
     etable &lt;- topTags(lrt, n=nrow(dge))\$table
     etable &lt;- etable[order(etable\$FDR), ]
@@ -175,7 +174,7 @@
     abline(h=c(-1,1), col="blue")
     dev.off()
   }
-  print("Done!")
+  write("Done!",stdout())
 }
 		</configfile>
 	</configfiles>