Mercurial > repos > yufei-luo > s_mart
view SMART/DiffExpAnal/wrappGSNAP.py @ 31:0ab839023fe4
Uploaded
author | m-zytnicki |
---|---|
date | Tue, 30 Apr 2013 14:33:21 -0400 |
parents | 94ab73e8a190 |
children |
line wrap: on
line source
#! /usr/bin/env python import os, sys, subprocess,tempfile from optparse import OptionParser def stop_err(msg): sys.stderr.write('%s\n' % msg) sys.exit() def __main__(): #Parse Command Line description = "GMAP/GSNAP version:2012-12-20." parser = OptionParser(description = description) parser.add_option("-d", "--genomeName", dest="genomeName", help="Define the reference genome name.[compulsory]") parser.add_option("-o", "--outputFile", dest="outputfile", help="output[compulsory]") #parser.add_option("-D", "--workingDir", dest="workingdir", help="Define the directory of writing reference genome index.[compulsory]") parser.add_option("-k", "--kmer", dest="kmer", default=12, help="Choose kmer value (<=16), a big kmer value can take more RAM(4Go).[compulsory]") parser.add_option("-i", "--inputFasta", dest="inputFastaFile", help="Reference genome file, fasta format.[compulsory]") parser.add_option("-q", "--inputFastq", dest="inputFastqFile", help="Input fastq file.") parser.add_option("-p", "--pairedEnd", dest="pairedEndFile", default=None, help="Input paired-end fastq file.") parser.add_option("-A", "--outputFormat", dest="outputFormat", default="sam", help="Choose an output format [sam, goby (need to re-compile with appropriate options)].") (options, args) = parser.parse_args() #If workingDir dose not exist, should create before run the job. workingDir = os.path.dirname(options.inputFastaFile) cmds = [] cmd_setup = "gmap_setup -d %s -D %s -k %s %s" % (options.genomeName, workingDir, options.kmer, options.inputFastaFile) cmds.append(cmd_setup) cmd_make_coords = "make -f Makefile.%s coords" % options.genomeName cmds.append(cmd_make_coords) cmd_make_gmapdb = "make -f Makefile.%s gmapdb" % options.genomeName cmds.append(cmd_make_gmapdb) cmd_make_install = "make -f Makefile.%s install" % options.genomeName cmds.append(cmd_make_install) cmd_run = "gsnap -d %s -D %s -A %s %s " % (options.genomeName, workingDir, options.outputFormat, options.inputFastqFile) if options.pairedEndFile != None: cmd_run += "%s" % options.pairedEndFile cmd_run += " > %s" % options.outputfile cmds.append(cmd_run) tmp_files = [] for i in range(len(cmds)): try: tmp_out = tempfile.NamedTemporaryFile().name tmp_files.append(tmp_out) tmp_stdout = open(tmp_out, 'wb') tmp_err = tempfile.NamedTemporaryFile().name tmp_files.append(tmp_err) tmp_stderr = open(tmp_err, 'wb') proc = subprocess.Popen(args=cmds[i], shell=True, cwd=".", stdout=tmp_stdout, stderr=tmp_stderr) returncode = proc.wait() tmp_stderr.close() #get stderr, allowing for case where it's very large tmp_stderr = open(tmp_err, 'rb') stderr = '' buffsize = 1048576 try: while True: stderr += tmp_stderr.read(buffsize) if not stderr or len(stderr) % buffsize != 0: break except OverflowError: pass tmp_stdout.close() tmp_stderr.close() if returncode != 0: raise Exception, stderr except Exception, e: stop_err('Error in :\n' + str(e)) for tmp_file in tmp_files: os.remove(tmp_file) if __name__=="__main__":__main__()