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