Mercurial > repos > fubar > digital_dge

# 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()
author	fubar
date	Fri, 09 Sep 2011 01:11:38 -0400
parents	1959becd0592
children