diff SMART/DiffExpAnal/wrappGSNAP.py @ 18:94ab73e8a190

Uploaded

--- /dev/null	Thu Jan 01 00:00:00 1970 +0000
+++ b/SMART/DiffExpAnal/wrappGSNAP.py	Mon Apr 29 03:20:15 2013 -0400
@@ -0,0 +1,76 @@
+#! /usr/bin/env python
+
+import os, sys, subprocess,tempfile
+from optparse import OptionParser
+
+def stop_err(msg):
+    sys.stderr.write('%s\n' % msg)
+    sys.exit()
+
+def __main__():
+    #Parse Command Line
+    description = "GMAP/GSNAP version:2012-12-20."
+    parser = OptionParser(description = description)
+    parser.add_option("-d", "--genomeName", dest="genomeName", help="Define the reference genome name.[compulsory]")
+    parser.add_option("-o", "--outputFile", dest="outputfile", help="output[compulsory]")
+    #parser.add_option("-D", "--workingDir", dest="workingdir", help="Define the directory of writing reference genome index.[compulsory]")
+    parser.add_option("-k", "--kmer", dest="kmer", default=12, help="Choose kmer value (<=16), a big kmer value can take more RAM(4Go).[compulsory]")
+    parser.add_option("-i", "--inputFasta", dest="inputFastaFile", help="Reference genome file, fasta format.[compulsory]")
+    parser.add_option("-q", "--inputFastq", dest="inputFastqFile", help="Input fastq file.")
+    parser.add_option("-p", "--pairedEnd", dest="pairedEndFile", default=None, help="Input paired-end fastq file.")
+    parser.add_option("-A", "--outputFormat", dest="outputFormat", default="sam", help="Choose an output format [sam, goby (need to re-compile with appropriate options)].")
+    (options, args) = parser.parse_args()    
+
+    #If workingDir dose not exist, should create before run the job.
+    
+    workingDir = os.path.dirname(options.inputFastaFile)
+    
+    cmds = []
+    cmd_setup = "gmap_setup -d %s -D %s -k %s %s" % (options.genomeName, workingDir, options.kmer, options.inputFastaFile)
+    cmds.append(cmd_setup)
+    cmd_make_coords = "make -f Makefile.%s coords" % options.genomeName 
+    cmds.append(cmd_make_coords)
+    cmd_make_gmapdb = "make -f Makefile.%s gmapdb" % options.genomeName
+    cmds.append(cmd_make_gmapdb)
+    cmd_make_install = "make -f Makefile.%s install" % options.genomeName
+    cmds.append(cmd_make_install)
+    cmd_run = "gsnap -d %s -D %s -A %s %s " % (options.genomeName, workingDir, options.outputFormat, options.inputFastqFile)
+    if options.pairedEndFile != None:
+        cmd_run += "%s" % options.pairedEndFile
+    cmd_run += " > %s" % options.outputfile
+    cmds.append(cmd_run)
+    
+    tmp_files = []
+    for i in range(len(cmds)):
+        try:
+            tmp_out = tempfile.NamedTemporaryFile().name
+            tmp_files.append(tmp_out)
+            tmp_stdout = open(tmp_out, 'wb')
+            tmp_err = tempfile.NamedTemporaryFile().name
+            tmp_files.append(tmp_err)
+            tmp_stderr = open(tmp_err, 'wb')
+            proc = subprocess.Popen(args=cmds[i], shell=True, cwd=".", stdout=tmp_stdout, stderr=tmp_stderr)
+            returncode = proc.wait()
+            tmp_stderr.close()
+            #get stderr, allowing for case where it's very large
+            tmp_stderr = open(tmp_err, 'rb')
+            stderr = ''
+            buffsize = 1048576
+            try:
+                while True:
+                    stderr += tmp_stderr.read(buffsize)
+                    if not stderr or len(stderr) % buffsize != 0:
+                        break
+            except OverflowError:
+                pass
+            tmp_stdout.close()
+            tmp_stderr.close()
+            if returncode != 0:
+                raise Exception, stderr
+        except Exception, e:
+            stop_err('Error in :\n' + str(e))
+    
+    for tmp_file in tmp_files:
+        os.remove(tmp_file)
+    
+if __name__=="__main__":__main__()