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