Mercurial > repos > yufei-luo > s_mart

comparison SMART/DiffExpAnal/wrappGSNAP.py @ 31:0ab839023fe4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author m-zytnicki
date Tue, 30 Apr 2013 14:33:21 -0400
parents 94ab73e8a190
children
comparison
equal deleted inserted replaced
30:5677346472b5 31:0ab839023fe4
1 #! /usr/bin/env python
2
3 import os, sys, subprocess,tempfile
4 from optparse import OptionParser
5
6 def stop_err(msg):
7 sys.stderr.write('%s\n' % msg)
8 sys.exit()
9
10 def __main__():
11 #Parse Command Line
12 description = "GMAP/GSNAP version:2012-12-20."
13 parser = OptionParser(description = description)
14 parser.add_option("-d", "--genomeName", dest="genomeName", help="Define the reference genome name.[compulsory]")
15 parser.add_option("-o", "--outputFile", dest="outputfile", help="output[compulsory]")
16 #parser.add_option("-D", "--workingDir", dest="workingdir", help="Define the directory of writing reference genome index.[compulsory]")
17 parser.add_option("-k", "--kmer", dest="kmer", default=12, help="Choose kmer value (<=16), a big kmer value can take more RAM(4Go).[compulsory]")
18 parser.add_option("-i", "--inputFasta", dest="inputFastaFile", help="Reference genome file, fasta format.[compulsory]")
19 parser.add_option("-q", "--inputFastq", dest="inputFastqFile", help="Input fastq file.")
20 parser.add_option("-p", "--pairedEnd", dest="pairedEndFile", default=None, help="Input paired-end fastq file.")
21 parser.add_option("-A", "--outputFormat", dest="outputFormat", default="sam", help="Choose an output format [sam, goby (need to re-compile with appropriate options)].")
22 (options, args) = parser.parse_args()
23
24 #If workingDir dose not exist, should create before run the job.
25
26 workingDir = os.path.dirname(options.inputFastaFile)
27
28 cmds = []
29 cmd_setup = "gmap_setup -d %s -D %s -k %s %s" % (options.genomeName, workingDir, options.kmer, options.inputFastaFile)
30 cmds.append(cmd_setup)
31 cmd_make_coords = "make -f Makefile.%s coords" % options.genomeName
32 cmds.append(cmd_make_coords)
33 cmd_make_gmapdb = "make -f Makefile.%s gmapdb" % options.genomeName
34 cmds.append(cmd_make_gmapdb)
35 cmd_make_install = "make -f Makefile.%s install" % options.genomeName
36 cmds.append(cmd_make_install)
37 cmd_run = "gsnap -d %s -D %s -A %s %s " % (options.genomeName, workingDir, options.outputFormat, options.inputFastqFile)
38 if options.pairedEndFile != None:
39 cmd_run += "%s" % options.pairedEndFile
40 cmd_run += " > %s" % options.outputfile
41 cmds.append(cmd_run)
42
43 tmp_files = []
44 for i in range(len(cmds)):
45 try:
46 tmp_out = tempfile.NamedTemporaryFile().name
47 tmp_files.append(tmp_out)
48 tmp_stdout = open(tmp_out, 'wb')
49 tmp_err = tempfile.NamedTemporaryFile().name
50 tmp_files.append(tmp_err)
51 tmp_stderr = open(tmp_err, 'wb')
52 proc = subprocess.Popen(args=cmds[i], shell=True, cwd=".", stdout=tmp_stdout, stderr=tmp_stderr)
53 returncode = proc.wait()
54 tmp_stderr.close()
55 #get stderr, allowing for case where it's very large
56 tmp_stderr = open(tmp_err, 'rb')
57 stderr = ''
58 buffsize = 1048576
59 try:
60 while True:
61 stderr += tmp_stderr.read(buffsize)
62 if not stderr or len(stderr) % buffsize != 0:
63 break
64 except OverflowError:
65 pass
66 tmp_stdout.close()
67 tmp_stderr.close()
68 if returncode != 0:
69 raise Exception, stderr
70 except Exception, e:
71 stop_err('Error in :\n' + str(e))
72
73 for tmp_file in tmp_files:
74 os.remove(tmp_file)
75
76 if __name__=="__main__":__main__()