comparison tools/ilmn_pacbio/quake_wrapper.py @ 0:9071e359b9a3

Uploaded
author xuebing
date Fri, 09 Mar 2012 19:37:19 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:9071e359b9a3
1 #!/usr/bin/python
2 #
3 # Copyright (c) 2011, Pacific Biosciences of California, Inc.
4 #
5 # All rights reserved.
6 #
7 #Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:
8 # * Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
9 # * Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
10 # * Neither the name of Pacific Biosciences nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.
11 #
12 #THIS SOFTWARE IS PROVIDED BY PACIFIC BIOSCIENCES AND ITS CONTRIBUTORS "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
13 #WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL PACIFIC BIOSCIENCES OR ITS CONTRIBUTORS BE LIABLE FOR ANY
14 #DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
15 #LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
16 #(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
17 #
18 import sys
19 import os
20 import subprocess
21
22 QUAKE_EXE = os.path.join( os.path.dirname(os.path.abspath(sys.argv[0])), 'quake.py' )
23 cmdLine = sys.argv
24 cmdLine.pop(0)
25
26 #
27 # horribly not robust, but it was a pain to rewrite everything with
28 # optparse
29 #
30 j = -1
31 cut = 0
32 for i,arg in enumerate(cmdLine):
33 if '--default_cutoff' in arg:
34 j = i
35 cut = int(arg.split('=')[1])
36 if j>=0:
37 cmdLine = cmdLine[:j] + cmdLine[j+1:]
38
39 j = -1
40 output=''
41 for i,arg in enumerate(cmdLine):
42 if '--output' in arg:
43 j = i
44 output = arg.split('=')[1]
45 if j>=0:
46 cmdLine = cmdLine[:j] + cmdLine[j+1:]
47
48 def backticks( cmd, merge_stderr=True ):
49 """
50 Simulates the perl backticks (``) command with error-handling support
51 Returns ( command output as sequence of strings, error code, error message )
52 """
53 if merge_stderr:
54 _stderr = subprocess.STDOUT
55 else:
56 _stderr = subprocess.PIPE
57
58 p = subprocess.Popen( cmd, shell=True, stdin=subprocess.PIPE,
59 stdout=subprocess.PIPE, stderr=_stderr,
60 close_fds=True )
61
62 out = [ l[:-1] for l in p.stdout.readlines() ]
63
64 p.stdout.close()
65 if not merge_stderr:
66 p.stderr.close()
67
68 # need to allow process to terminate
69 p.wait()
70
71 errCode = p.returncode and p.returncode or 0
72 if p.returncode>0:
73 errorMessage = os.linesep.join(out)
74 output = []
75 else:
76 errorMessage = ''
77 output = out
78
79 return output, errCode, errorMessage
80
81 def to_stdout():
82 def toCorFastq(f):
83 stem, ext = os.path.splitext( os.path.basename(f) )
84 dir = os.path.dirname(f)
85 corFastq = os.path.join(dir,'%s.cor%s' % (stem,ext) )
86 if not os.path.exists(corFastq):
87 print >>sys.stderr, "Can't find path %s" % corFastq
88 sys.exit(1)
89 return corFastq
90 if '-r' in cmdLine:
91 fastqFile = cmdLine[ cmdLine.index('-r')+1 ]
92 corFastq = toCorFastq(fastqFile)
93 infile = open( corFastq, 'r' )
94 for line in infile:
95 sys.stdout.write( line )
96 infile.close()
97 else:
98 fofnFile = cmdLine[ cmdLine.index('-f')+1 ]
99 infile = open(fofnFile,'r')
100 for line in infile:
101 line = line.strip()
102 if len(line)>0:
103 fastqFiles = line.split()
104 break
105 infile.close()
106 outs = output.split(',')
107 for o,f in zip(outs,fastqFiles):
108 cf = toCorFastq(f)
109 os.system( 'cp %s %s' % ( cf, o ) )
110
111 def run():
112 cmd = '%s %s' % ( QUAKE_EXE, " ".join(cmdLine) )
113 output, errCode, errMsg = backticks( cmd )
114
115 if errCode==0:
116 to_stdout()
117 else:
118 # if Quake exits with an error in cutoff determination we
119 # can force correction if requested
120 if 'cutoff.txt' in errMsg and cut>0:
121 outfile = open( 'cutoff.txt', 'w' )
122 print >>outfile, str(cut)
123 outfile.close()
124 cmd = '%s --no_count --no_cut %s' % ( QUAKE_EXE, " ".join(cmdLine) )
125 output, errCode, errMsg = backticks( cmd )
126 if errCode==0:
127 to_stdout()
128 else:
129 print >>sys.stderr, errMsg
130 sys.exit(1)
131
132 if __name__=='__main__': run()