|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 import os, sys, tempfile
|
|
|
4
|
|
|
5 assert sys.version_info[:2] >= (2.4)
|
|
|
6
|
|
|
7 def stop_err( msg ):
|
|
|
8
|
|
|
9 sys.stderr.write( "%s\n" % msg )
|
|
|
10 sys.exit()
|
|
|
11
|
|
|
12
|
|
|
13 def __main__():
|
|
|
14
|
|
|
15 # I/O
|
|
|
16 target_path = sys.argv[1]
|
|
|
17 infile = sys.argv[2]
|
|
|
18 scorefile = sys.argv[3]
|
|
|
19 high_score = sys.argv[4] # -q
|
|
|
20 high_len = sys.argv[5] # -M
|
|
|
21 read_len = sys.argv[6] # -w
|
|
|
22 align_len = sys.argv[7] # -h
|
|
|
23 mismatch = sys.argv[8] # -m
|
|
|
24 output_file = sys.argv[9]
|
|
|
25
|
|
|
26 try:
|
|
|
27 float(high_score)
|
|
|
28 except:
|
|
|
29 stop_err('Invalid value for minimal quality score.')
|
|
|
30
|
|
|
31 try:
|
|
|
32 int(high_len)
|
|
|
33 except:
|
|
|
34 stop_err('Invalid value for minimal high quality bases.')
|
|
|
35
|
|
|
36 # first guess the read length
|
|
|
37 guess_read_len = 0
|
|
|
38 seq = ''
|
|
|
39 for i, line in enumerate(open(infile)):
|
|
|
40 line = line.rstrip('\r\n')
|
|
|
41 if line.startswith('>'):
|
|
|
42 if seq:
|
|
|
43 guess_read_len = len(seq)
|
|
|
44 break
|
|
|
45 else:
|
|
|
46 seq += line
|
|
|
47
|
|
|
48 try:
|
|
|
49 test = int(read_len)
|
|
|
50 if test == 0:
|
|
|
51 read_len = str(guess_read_len)
|
|
|
52 else:
|
|
|
53 assert test >= 20 and test <= 64
|
|
|
54 except:
|
|
|
55 stop_err('Invalid value for read length. Must be between 20 and 64.')
|
|
|
56
|
|
|
57
|
|
|
58 try:
|
|
|
59 int(align_len)
|
|
|
60 except:
|
|
|
61 stop_err('Invalid value for minimal length of a hit.')
|
|
|
62
|
|
|
63 try:
|
|
|
64 int(mismatch)
|
|
|
65 except:
|
|
|
66 stop_err('Invalid value for mismatch numbers in an alignment.')
|
|
|
67
|
|
|
68 all_files = []
|
|
|
69 if os.path.isdir(target_path):
|
|
|
70 # check target genome
|
|
|
71 fa_files = os.listdir(target_path)
|
|
|
72
|
|
|
73 for file in fa_files:
|
|
|
74 file = "%s/%s" % ( target_path, file )
|
|
|
75 file = os.path.normpath(file)
|
|
|
76 all_files.append(file)
|
|
|
77 else:
|
|
|
78 stop_err("No sequences for %s are available for search, please report this error." %(target_path))
|
|
|
79
|
|
|
80 for detail_file_path in all_files:
|
|
|
81 output_tempfile = tempfile.NamedTemporaryFile().name
|
|
|
82 command = "rmapq -q %s -M %s -h %s -w %s -m %s -Q %s -c %s %s -o %s 2>&1" % ( high_score, high_len, align_len, read_len, mismatch, scorefile, detail_file_path, infile, output_tempfile )
|
|
|
83 #print command
|
|
|
84 try:
|
|
|
85 os.system( command )
|
|
|
86 except Exception, e:
|
|
|
87 stop_err( str( e ) )
|
|
|
88
|
|
|
89 try:
|
|
|
90 assert os.system( 'cat %s >> %s' % ( output_tempfile, output_file ) ) == 0
|
|
|
91 except Exception, e:
|
|
|
92 stop_err( str( e ) )
|
|
|
93
|
|
|
94 try:
|
|
|
95 os.remove( output_tempfile )
|
|
|
96 except:
|
|
|
97 pass
|
|
|
98
|
|
|
99
|
|
|
100 if __name__ == '__main__': __main__()
|
