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

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