10
|
1 #!/usr/bin/python
|
|
2
|
|
3 import argparse, optparse, os, shutil, subprocess, sys, tempfile, shlex, time
|
|
4
|
|
5
|
|
6 parser = argparse.ArgumentParser(description='')
|
|
7
|
|
8 parser.add_argument ( '-i1', dest='input_bam', help='the bam input file' )
|
|
9 parser.add_argument ( '-o1', dest='output_raw', help='the output file' )
|
|
10 parser.add_argument ( '-o2', dest='output_vcf', help='the output file' )
|
|
11
|
|
12
|
|
13 tmp_dir = tempfile.mkdtemp()
|
|
14 errorFile = tmp_dir+"/errorLog"
|
|
15
|
16
|
16 currentDir = os.path.dirname(os.path.abspath(__file__))
|
21
|
17
|
22
|
18 binPath = os.environ['BREAKDANCER_BIN']
|
25
|
19 print(binPath)
|
|
20 print(str(binPath))
|
|
21 bam2cfg = str(binPath)+"/bam2cfg.pl"
|
21
|
22 breakdancer = binPath+"/breakdancer-max"
|
|
23 breakdancer2vcf = binPath+"/breakdancer2vcf.py"
|
10
|
24
|
|
25
|
|
26 def bam2cfg(args):
|
|
27 config = tmp_dir+"/breakdancer_config"
|
19
|
28 cmd = 'perl %s %s' % (bam2cfg, args.input_bam)
|
10
|
29 execute(cmd, output=config)
|
|
30 print ("\ncmd = %s \n" ) %(cmd)
|
|
31 return config
|
|
32
|
|
33
|
|
34 def breakdancer(args, config):
|
19
|
35 cmd = '%s %s' % (breakdancer, config)
|
10
|
36 execute(cmd, output=args.output_raw)
|
|
37
|
|
38
|
|
39 def breakdancer2vcf(args):
|
19
|
40 cmd = "python %s -i %s -o %s" % ( breakdancer2vcf, args.output_raw, args.output_vcf )
|
10
|
41 execute(cmd)
|
|
42
|
|
43
|
|
44 def execute( cmd, output="" ):
|
|
45 try:
|
|
46 err = open(tmp_dir+"/errorLog", 'a')
|
|
47 if output != "":
|
|
48 out = open(output, 'w')
|
|
49 else:
|
|
50 out = subprocess.PIPE
|
|
51 process = subprocess.Popen( args=shlex.split(cmd), stdout=out, stderr=err )
|
|
52 process.wait()
|
|
53 err.close()
|
|
54 if out != subprocess.PIPE:
|
|
55 out.close()
|
|
56 except Exception, e:
|
|
57 sys.stderr.write("problem doing : %s\n" %(cmd))
|
|
58 sys.stderr.write( '%s\n\n' % str(e) )
|
|
59
|
|
60
|
|
61 def check( output ):
|
|
62 if (os.path.getsize(output)>0):
|
|
63 return True
|
|
64 else:
|
|
65 sys.stderr.write('The output file is empty : %s\n' % (output))
|
|
66
|
|
67
|
|
68 def getLine(file):
|
|
69 try:
|
|
70 f = open(file, 'r')
|
|
71 lignes = f.readlines()
|
|
72 n=0
|
|
73 for ligne in lignes:
|
|
74 if ligne.strip()[0]!="#":
|
|
75 n+=1
|
|
76 # n = len(lignes)
|
|
77 f.close()
|
|
78 except Exception, e:
|
|
79 sys.stderr.write( '%s\n' % str(e) )
|
|
80 return n
|
|
81
|
|
82
|
|
83 def compare( output1, output2 ):
|
|
84 # compare le nombre de ligne entre 2 fichiers
|
|
85 # pour verifier qu'aucune ligne n'a ete saute
|
|
86 num_raw = getLine(output1)
|
|
87 num_vcf = getLine(output2)
|
|
88 if (num_raw==num_vcf):
|
|
89 return True
|
|
90 else:
|
|
91 sys.stderr.write('Not the same number of variant between the raw file and the vcf file : %d vs %d\n' % (num_raw, num_vcf))
|
|
92
|
|
93
|
|
94 def __main__():
|
|
95
|
|
96 time0 = time.time()
|
|
97 args = parser.parse_args()
|
|
98
|
|
99 try:
|
|
100 config = bam2cfg(args)
|
|
101 breakdancer(args, config)
|
|
102 breakdancer2vcf(args)
|
|
103
|
|
104 # quelques tests
|
|
105 # verifier que les fichiers de sorties ne sont pas vides
|
|
106 check(args.output_raw)
|
|
107 check(args.output_vcf)
|
|
108 # comparer le nombre de ligne entre les 2 fichiers
|
|
109 # pour etre sur que toute les variations ont ete prises en compte
|
|
110 compare(args.output_raw, args.output_vcf)
|
|
111
|
|
112 sys.stdout.write( '\nDone in %d seconds\n' %(int(time.time()-time0)))
|
|
113
|
|
114 if ( os.path.getsize(errorFile)>0 ):
|
|
115 sys.stdout.write( 'At least one non fatal error was send, check the file : %s\n' %(errorFile) )
|
|
116 try:
|
|
117 err = open(errorFile, 'r')
|
|
118 errors = err.read()
|
|
119 err.close()
|
|
120 sys.stdout.write( "Errors :\n%s" %(errors))
|
|
121 except Exception, e:
|
|
122 sys.stderr.write( '%s\n' % str(e) )
|
|
123 else:
|
|
124 sys.stdout.write( 'BreakDancer successful' )
|
|
125
|
|
126 except Exception, e:
|
|
127 sys.stdout.write( 'BreakDancer fail\n' )
|
|
128 sys.stderr.write( '%s\n' % str(e) )
|
|
129 sys.exit()
|
|
130
|
|
131 finally:
|
|
132 if os.path.exists( tmp_dir ):
|
|
133 shutil.rmtree( tmp_dir )
|
|
134
|
|
135 if __name__=="__main__":
|
|
136 __main__()
|