Mercurial > repos > fubar > digital_dge

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/rgDGE.py	Fri Sep 09 01:07:48 2011 -0400
@@ -0,0 +1,338 @@
+# Copyright ross lazarus
+# september 2011
+# all rights reserved
+# for the Rgenetics project
+# all rights reserved
+# licensed to you under the terms of the LGPL as documented at http://www.gnu.org/copyleft/lesser.html
+
+import sys
+import shutil
+import subprocess
+import os
+import time
+import tempfile
+import optparse
+
+progname = os.path.split(sys.argv[0])[1]
+myversion = 'V000.2 September 2011'
+verbose = False
+debug = False
+
+
+galhtmlprefix = """<?xml version="1.0" encoding="utf-8" ?>
+<!DOCTYPE html PUBLIC "-//W3C//DTD XHTML 1.0 Transitional//EN" "http://www.w3.org/TR/xhtml1/DTD/xhtml1-transitional.dtd">
+<html xmlns="http://www.w3.org/1999/xhtml" xml:lang="en" lang="en">
+<head> <meta http-equiv="Content-Type" content="text/html; charset=utf-8" />
+<meta name="generator" content="Galaxy %s tool output - see http://g2.trac.bx.psu.edu/" />
+<title></title>
+<link rel="stylesheet" href="/static/style/base.css" type="text/css" />
+</head>
+<body>
+<div class="document">
+"""
+galhtmlattr = """<b><a href="http://rgenetics.org">Galaxy Rgenetics</a> tool output %s run at %s</b><br/>"""
+galhtmlpostfix = """</div></body></html>\n"""
+
+DGESCRIPT="""
+#### edgeR.Rscript
+#
+# Performs DGE on a count table containing n replicates of two conditions
+#
+# Parameters
+#
+# 1 - Output Dir
+
+# Writen by: S.Lunke and A.Kaspi
+####
+
+
+usage <- function(){
+  print("#### edgeR.R", quote=F)
+  print("", quote=F)
+  print("Performs DGE on a count table containing n replicates of two conditions", quote=F)
+  print("", quote=F)
+  print("USAGE: Rscript edgeR.R <OUT_DIR> <INPUT> <TreatmentName> <ControlName> <Treatcol1,2,3> <Controlcol1,2,3>", quote=F)
+  print("", quote=F)
+  print(" Parameters", quote=F)
+  print("", quote=F)
+  print(" 1 - Output Dir", quote=F)
+  print(" 2 - Input File", quote=F)
+  print(" 3 - Treatment name", quote=F)
+  print(" 4 - Treatment Columns i.e. 3,4,6", quote=F)
+  print(" 5 - Control name", quote=F)
+  print(" 6 - Control Columns i.e. 2,7,8", quote=F)
+  print(" 7 - Output tabular file name", quote=F)
+
+
+  print("", quote=F)
+  print(" Interface writen by: S.Lunke and A.Kaspi", quote=F)
+  print(" Makes extensive use of the edgeR software - Mark D. Robinson, Davis J. McCarthy, Gordon K. Smyth, PMCID: PMC2796818", quote=F)
+  print("####", quote=F)
+  q()
+   }
+
+
+## edgeIt runs edgeR
+edgeIt <- function (Input,group,outputfilename) {
+
+
+        ## Error handling
+        if (length(unique(group))!=2){
+                print("Number of conditions identified in experiment does not equal 2")
+                q()
+        }
+
+        #The workhorse
+        require(edgeR)
+
+        ## Setup DGEList object
+        DGEList <- DGEList(Count_Matrix, group = group)
+
+        #Extract T v C names
+        TName=unique(group)[1]
+        CName=unique(group)[2]
+
+        # Outfile name - we need something more predictable than
+        # outfile <- paste(Out_Dir,"/",TName,"_Vs_",CName, sep="")
+        # since it needs to be renamed and squirted into the history so added as a paramter
+
+        # Filter all tags that have less than one read per million in half + 1 or more libraries. n = ceiling(length(colnames(DGEList))/2)+1
+        n = ceiling(length(colnames(DGEList))/2)+1
+        DGEList <- DGEList[rowSums(1e+06 * DGEList$counts/expandAsMatrix(DGEList$samples$lib.size, dim(DGEList)) > 1) >= n, ]
+
+
+        ## Run tagwise dispersion test and DGE analysis
+        prior.n <- ceiling(50/(length(group)-length(2)))
+        DGEList <- calcNormFactors(DGEList)
+        DGEList <- estimateCommonDisp(DGEList)
+        DGEList <- estimateTagwiseDisp(DGEList, prior.n=prior.n, trend=T, grid.length=1000)
+        DE <- exactTest(DGEList, common.disp=F)
+
+        ## Normalized data RPM
+        normData <- (1e+06 * DGEList$counts/expandAsMatrix(DGEList$samples$lib.size, dim(DGEList)))
+        colnames(normData) <- paste( "norm",colnames(normData),sep="_")
+
+
+        #Prepare our output file
+        output <- cbind(
+                Name=as.character(rownames(DGEList$counts)),
+                DE$table,
+                adj.p.value=p.adjust(DE$table$p.value, method="fdr"),
+                Dispersion=DGEList$tagwise.dispersion,normData,
+                DGEList$counts
+        )
+        soutput = output[order(output$p.val),] # sorted into p value order - for quick toptable
+        #Write output
+        write.table(soutput,outputfilename, quote=FALSE, sep="\t",row.names=F)
+
+
+        ## Plot MAplot
+        pdf("Smearplot.pdf")
+        #TODO colour siggenes
+        nsig = sum(output$adj.p.value < 0.05)
+        deTags = rownames(topTags(DE,n=nsig)$table)
+        plotSmear(DGEList,de.tags=deTags,main=paste("Smear Plot for ",TName,' Vs ',CName,' (FDR@0.05, N = ',nsig,')',sep=''))
+        grid(col="blue")
+        dev.off()
+
+        ## Plot MDS
+        pdf("MDSplot.pdf")
+        plotMDS.dge(DGEList,xlim=c(-2,1),main=paste("MDS Plot for",TName,'Vs',CName),cex=0.5)
+        grid(col="blue")
+        dev.off()
+
+        #Return our main table
+        output
+
+}       #Done
+
+#### MAIN ####
+
+parameters <- commandArgs(trailingOnly=T)
+
+## Error handling
+if (length(parameters) < 6){
+  print("Not enough Input files supplied. Specify at least two input files.", quote=F)
+  print("", quote=F)
+  usage()
+}
+
+Out_Dir         <- as.character(parameters[1])
+Input           <- as.character(parameters[2])
+TreatmentName   <- as.character(parameters[3])
+TreatmentCols   <- as.character(parameters[4])
+ControlName     <- as.character(parameters[5])
+ControlCols     <- as.character(parameters[6])
+outputfilename  <- as.character(parameters[7])
+
+
+#Set our columns
+TCols           = as.numeric(strsplit(TreatmentCols,",")[[1]])-1 #+1
+CCols           = as.numeric(strsplit(ControlCols,",")[[1]])-1 #+1
+cat('## got TCols=')
+cat(TCols)
+cat('; CCols=')
+cat(CCols)
+cat('\n')
+
+
+group<-strsplit(TreatmentCols,",")[[1]]
+
+## Create output dir if non existent
+  if (file.exists(Out_Dir) == F) dir.create(Out_Dir)
+
+Count_Matrix<-read.table(Input,header=T,row.names=1,sep='\t')                           #Load tab file assume header
+Count_Matrix<-Count_Matrix[,c(TCols,CCols)]                                             #extract columns we care about
+group<-c(rep(TreatmentName,length(TCols)), rep(ControlName,length(CCols)) )             #Build a group descriptor
+colnames(Count_Matrix) <- paste(group,colnames(Count_Matrix),sep="_")                   #Relable columns
+
+
+results <- edgeIt(Input,group,outputfilename)                                           #Run the main function
+# for the log
+sessionInfo()
+
+"""
+
+def timenow():
+    """return current time as a string
+    """
+    return time.strftime('%d/%m/%Y %H:%M:%S', time.localtime(time.time()))
+
+class DGE:
+    """class is a wrapper for DGE - note hard coded script above so it's all in one place for Galaxy"""
+
+    def __init__(self,myName=None,opts=None):
+        """
+        Rscript edgeR.R results HGHN_mRNA_counts.txt HiGly 2,3 Control 1,4
+Out_Dir         <- as.character(parameters[1])
+Input           <- as.character(parameters[2])
+TreatmentName   <- as.character(parameters[3])
+TreatmentCols   <- as.character(parameters[4])
+ControlName     <- as.character(parameters[5])
+ControlCols     <- as.character(parameters[6])
+outputfilename  <- as.character(parameters[7])
+        """
+        self.thumbformat = 'jpg'
+        self.tlog = os.path.join(opts.output_dir,"DGE_runner.log")
+        self.opts = opts
+        self.myName=myName
+        self.cl = []
+        a = self.cl.append
+        rfname = os.path.join(opts.output_dir,'DGE.R')
+        rscript = open(rfname,'w')
+        rscript.write(DGESCRIPT)
+        rscript.close()
+        a('Rscript')
+        a(rfname)
+        a("%s" % opts.output_dir)
+        a("%s" % opts.input_matrix)
+        a("%s" % opts.treatname)
+        a("%s" % opts.treatcols)
+        a("%s" % opts.ctrlname)
+        a("%s" % opts.ctrlcols)
+        a("%s" % opts.output_tab)
+
+    def compressPDF(self,inpdf=None,thumbformat='png'):
+        """need absolute path to pdf
+        """
+        assert os.path.isfile(inpdf), "## Input %s supplied to %s compressPDF not found" % (inpdf,self.myName)
+        sto = open(self.tlog,'a')
+        outpdf = '%s_compressed' % inpdf
+        cl = ["gs", "-sDEVICE=pdfwrite", "-dNOPAUSE", "-dBATCH", "-sOutputFile=%s" % outpdf,inpdf]
+        x = subprocess.Popen(cl,stdout=sto,stderr=sto,cwd=self.opts.output_dir)
+        retval1 = x.wait()
+        if retval1 == 0:
+            os.unlink(inpdf)
+            shutil.move(outpdf,inpdf)
+        outpng = '%s.%s' % (os.path.splitext(inpdf)[0],thumbformat)
+        cl2 = ['convert', inpdf, outpng]
+        x = subprocess.Popen(cl2,stdout=sto,stderr=sto,cwd=self.opts.output_dir)
+        retval2 = x.wait()
+        sto.close()
+        retval = retval1 or retval2
+        return retval
+
+    def runDGE(self):
+        """
+        """
+        sto = open(self.tlog,'w')
+        x = subprocess.Popen(' '.join(self.cl),shell=True,stdout=sto,stderr=sto,cwd=self.opts.output_dir)
+        retval = x.wait()
+        sto.close()
+        flist = os.listdir(self.opts.output_dir)
+        flist.sort()
+        html = [galhtmlprefix % progname,]
+        html.append('<h2>Galaxy DGE outputs run at %s<h2></br>Click on the images below to download high quality PDF versions</br>\n' % timenow())
+        if len(flist) > 0:
+            html.append('<table>\n')
+            for fname in flist:
+                dname,e = os.path.splitext(fname)
+                if e.lower() == '.pdf' : # compress and make a thumbnail
+                    thumb = '%s.%s' % (dname,self.thumbformat)
+                    pdff = os.path.join(self.opts.output_dir,fname)
+                    retval = self.compressPDF(inpdf=pdff,thumbformat=self.thumbformat)
+                    if retval == 0:
+                        s= '<tr><td><a href="%s"><img src="%s" title="Click to download a print quality PDF of %s" hspace="10" width="600"></a></td></tr>' \
+                         % (fname,thumb,fname)
+                    else:
+                        dname = '%s (thumbnail image not_found)' % fname
+                        s= '<tr><td><a href="%s">%s</a></td></tr>' % (fname,dname)
+                    html.append(s)
+                else:
+                   html.append('<tr><td><a href="%s">%s</a></td></tr>' % (fname,fname))
+            html.append('</table>\n')
+        else:
+            html.append('<h2>### Error - R returned no files - please confirm that parameters are sane</h1>')
+        html.append('<h3>R log follows below</h3><hr><pre>\n')
+        rlog = open(self.tlog,'r').readlines()
+        html += rlog
+        html.append('%s CL = %s</br>\n' % (self.myName,' '.join(sys.argv)))
+        html.append('DGE.R CL = %s</br>\n' % (' '.join(self.cl)))
+        html.append('</pre>\n')
+        html.append(galhtmlattr % (progname,timenow()))
+        html.append(galhtmlpostfix)
+        htmlf = file(self.opts.output_html,'w')
+        htmlf.write('\n'.join(html))
+        htmlf.write('\n')
+        htmlf.close()
+
+
+def main():
+    u = """
+    This is a Galaxy wrapper. It expects to be called by DGE.xml as:
+    <command interpreter="python">rgDGE.py --output_html "$html_file.files_path" --input_matrix "$input1" --treatcols "$Treat_cols" --treatname "$treatment_name"
+--ctrlcols "$Control_cols"
+      --ctrlname "$control_name" --output_tab "$outtab" --output_html "$html_file" --output_dir "$html_file.files_path" --method "edgeR"
+    </command>
+    """
+    op = optparse.OptionParser()
+    a = op.add_option
+    a('--input_matrix',default=None)
+    a('--output_tab',default=None)
+    a('--output_html',default=None)
+    a('--treatcols',default=None)
+    a('--treatname',default='Treatment')
+    a('--ctrlcols',default=None)
+    a('--ctrlname',default='Control')
+    a('--output_dir',default=None)
+    a('--method',default='edgeR')
+    opts, args = op.parse_args()
+    assert opts.input_matrix and opts.output_html,u
+    assert os.path.isfile(opts.input_matrix),'## DGE runner unable to open supplied input file %s' % opts.input_matrix
+    assert opts.treatcols,'## DGE runner requires a comma separated list of treatment columns eg --treatcols 4,5,6'
+    assert opts.treatcols,'## DGE runner requires a comma separated list of control columns eg --ctlcols 2,3,7'
+    myName=os.path.split(sys.argv[0])[-1]
+    if not os.path.exists(opts.output_dir):
+        os.makedirs(opts.output_dir)
+    m = DGE(myName, opts=opts)
+    retcode = m.runDGE()
+    if retcode:
+        sys.exit(retcode) # indicate failure to job runner
+
+
+
+if __name__ == "__main__":
+    main()
+
+
+