| 18 | 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__() |