Mercurial > repos > jeremie > delly_

comparison delly.py @ 9:bccaed23a22b draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author jeremie
date Tue, 17 Jun 2014 08:27:47 -0400
parents
children
comparison
equal deleted inserted replaced
8:ab3c1999f607 9:bccaed23a22b
1 #!/usr/bin/python
2
3 import argparse, os, shutil, subprocess, sys, tempfile, shlex, vcf
4
5 parser = argparse.ArgumentParser(description='')
6 parser.add_argument ( '-b', dest='bam', help='the bam file', required=True )
7 parser.add_argument ( '-t', dest='types', help='SV analysis type (DEL, DUP, INV, TRA)', nargs='+', required=True )
8 parser.add_argument ( '-q', dest='map_qual', help='min. paired-end mapping quality', default='0' )
9 parser.add_argument ( '-s', dest='mad_cutoff', help='insert size cutoff, median+s*MAD (deletions only)', default=5 )
10 parser.add_argument ( '-g', dest='genome', help='genome fasta file' )
11 parser.add_argument ( '-m', dest='min_flank', help='minimum flanking sequence size', default=13 )
12 parser.add_argument ( '-e', dest='epsilon', help='allowed epsilon deviation of PE vs. SR deletion', default=0.1 )
13 parser.add_argument ( '-o', dest='output', help='output file' )
14
15 dellyPath="delly_v0.5.5"
16
17 def execute( cmd, output="" ):
18 tmp_dir = tempfile.mkdtemp()
19 try:
20 err = open(tmp_dir+"/errorLog", 'a')
21 if output != "":
22 out = open(output, 'w')
23 else:
24 out = subprocess.PIPE
25 process = subprocess.Popen( args=shlex.split(cmd), stdout=out, stderr=err )
26 process.wait()
27 err.close()
28 if out != subprocess.PIPE:
29 out.close()
30 except Exception, e:
31 sys.stderr.write("problem doing : %s\n" %(cmd))
32 sys.stderr.write( '%s\n\n' % str(e) )
33
34
35 def delly(args, tempDir):
36 tempOutputs=[];
37 for typ in args.types:
38 output=str(tempDir)+"/"+str(typ)
39 tempOutputs.append(output)
40 cmd = "%s %s -t %s -o %s -q %s -s %s -m %s -e %s " % (dellyPath, args.bam, typ, output, args.map_qual, args.mad_cutoff, args.min_flank, args.epsilon)
41 if (args.genome!="" and args.genome!=None):
42 cmd += "-g %s" % (args.genome)
43 # print ("commande = "+cmd)
44 execute( cmd )
45 return tempOutputs
46
47
48 def merge(outputs, outputFile):
49 template = vcf.Reader(filename=outputs[0])
50 vcfWriter = vcf.Writer(open(outputFile, 'w'), template)
51 for output in outputs:
52 vcfReader = vcf.Reader(filename=output)
53 for record in vcfReader:
54 vcfWriter.write_record(record)
55 return 0
56
57
58 def getVersion(program):
59 try:
60 tmp = tempfile.NamedTemporaryFile().name
61 tmp_stdout = open( tmp, 'wb' )
62 proc = subprocess.Popen( args=program, shell=True, stdout=tmp_stdout )
63 tmp_stdout.close()
64 returncode = proc.wait()
65 stdout = None
66 for line in open( tmp_stdout.name, 'rb' ):
67 if line.lower().find( 'version' ) >= 0:
68 stdout = line.strip()
69 break
70 if stdout:
71 sys.stdout.write( '%s\n' % stdout )
72 except:
73 sys.stdout.write( 'Could not determine %s version\n' % (program) )
74
75
76 def __main__():
77
78 args = parser.parse_args()
79
80 tempDir = tempfile.mkdtemp();
81 getVersion(dellyPath)
82
83 try:
84 execute("samtools index " + args.bam)
85 except Exception, e:
86 print("problem while indexing bam file " + str(e))
87
88 try:
89 tempOutputs = delly(args, tempDir)
90 except Exception, e:
91 sys.stdout.write("problem while runing delly " + str(e))
92
93 try:
94 if args.output:
95 outputs=[]
96 for output in tempOutputs:
97 if os.path.exists(output):
98 if os.path.getsize(output)>0:
99 outputs.append(output)
100 if len(outputs)>0:
101 merge(outputs, args.output)
102 except Exception, e:
103 sys.stdout.write("problem while merging files " + str(e))
104
105 finally:
106 try:
107 if os.path.exists(tempDir):
108 shutil.rmtree(tempDir)
109 os.system("rm "+str(args.bam)+".bai")
110 except:
111 pass
112
113 if __name__=="__main__":
114 __main__()