|
2
|
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
|
