Mercurial > repos > subazini > ngsaligners

comparison stampy_wrapper.py @ 2:55b3da30e4e5 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author subazini
date Wed, 17 Dec 2014 10:16:35 -0500
parents
children
comparison
equal deleted inserted replaced
1:abe73a62b59a 2:55b3da30e4e5
1 #!/usr/bin/env python
2
3 import optparse, os, shutil, subprocess, sys, tempfile, fileinput, commands
4 import time
5 import re
6 import psutil
7 def stop_err( msg ):
8 sys.stderr.write( "%s\n" % msg )
9 sys.exit()
10
11 def __main__():
12 #Parse Command Line
13 start = time.time()
14 #print "start time", start
15 parser = optparse.OptionParser()
16 parser.add_option( '-s','--species',dest='species')
17 parser.add_option( '-a', '--assembly', dest='assembly' )
18 parser.add_option( '-G', '--genome_index', dest='genome_index' )
19 parser.add_option( '-e', '--genome', dest='genome')
20 parser.add_option( '-D', '--indexes-path', dest='index_path', help='Indexes directory; location of .stidx and .sthash files.' )
21
22 # Wrapper options.
23 parser.add_option( '-O', '--output', dest='output' )
24 parser.add_option( '-1', '--input1', dest='input1', help='The (forward or single-end) reads file in Sanger FASTQ format' )
25 parser.add_option( '-2', '--input2', dest='input2', help='The reverse reads file in Sanger FASTQ format' )
26 parser.add_option( '-g', '--use_genome', dest='use_genome' )
27 parser.add_option( '-H', '--hashbuild', dest='build_hash' )
28 parser.add_option( '-M', '--map', dest='map', action="store_true" )
29
30 # Read group options.
31 parser.add_option( '-d', '--sd', dest='sd' )
32 parser.add_option( '-i', '--insert', dest='insert' )
33 parser.add_option( '-r', '--subrate', dest='subrate' )
34
35 parser.add_option( '-t', '--settings', dest='settings' )
36
37
38 (options, args) = parser.parse_args()
39
40 # Creat genome index if necessary.
41 tmp_index_dir = tempfile.mkdtemp()
42
43 # Genome index from history
44 if options.genome_index:
45 index_path = os.path.join( tmp_index_dir, '.'.join( os.path.split( options.genome_index )[1].split( '.' )[:-1] ) )
46 try:
47 os.link( options.genome_index, index_path + '.fa' )
48 #print options.genome_index
49 except:
50 pass
51 # Indexing genome
52 cmd_index = 'stampy.py --species=%s --assembly=%s -G %s %s' %(options.species,options.assembly,options.genome_index,options.genome)
53 #print cmd_index
54 #print options.genome
55
56 try:
57 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
58 tmp_stderr = open( tmp, 'wb' )
59 proc=subprocess.Popen(cmd_index,shell=True)
60 returncode = proc.wait()
61 tmp_stderr.close()
62 tmp_stderr = open( tmp, 'rb' )
63 stderr = ''
64 buffsize = 1048576
65 try:
66 while True:
67 stderr += tmp_stderr.read( buffsize )
68 if not stderr or len( stderr ) % buffsize != 0:
69 break
70 except OverflowError:
71 pass
72 tmp_stderr.close()
73
74 except Exception, e:
75 if os.path.exists( tmp_index_dir ):
76 shutil.rmtree( tmp_index_dir )
77 stop_err( 'Error indexing reference sequence\n' + str( e ) )
78 else:
79 index_path = options.index_path
80 try:
81 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir).name
82 tmp_stderr = open( tmp, 'wb' )
83 # Hashing genome
84 cmd_index2 = 'stampy.py -g %s -H %s > %s' %(options.genome_index,options.genome_index,options.output)
85 #print cmd_index2
86 proc=subprocess.Popen(cmd_index2, shell=True)
87 returncode = proc.wait()
88 #print proc.returncode
89
90 except Exception, e:
91 if os.path.exists( tmp_index_dir ):
92 shutil.rmtree( tmp_index_dir )
93 cmd_index3 = 'stampy.py -g %s -h %s' %(options.genome_index,options.genome_index)
94 if options.settings == "full":
95 opts = ' --insertsize=%s --insertsd=%s --substitutionrate=%s -M %s %s > %s' %(options.insert,options.sd,options.subrate,options.input1,options.input2,options.output)
96 else:
97 opts = '-M %s %s > %s' %(options.input1,options.input2,options.output)
98 cmd_index3 += opts
99 try:
100 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
101 tmp_stderr = open( tmp, 'wb' )
102 #print cmd_index3
103 proc=subprocess.call(cmd_index3, shell=True)
104 except Exception, e:
105 if os.path.exists( tmp_index_dir ):
106 shutil.rmtree( tmp_index_dir )
107
108 end = time.time()
109 #print "end time", end
110 elapsed = end - start
111 #print ("Time elapsed: {0:4f}sec",format(elapsed))
112
113 if __name__=="__main__": __main__()
114
115