0
|
1 #!/usr/bin/env python
|
|
2
|
|
3 """
|
|
4 Creates a pileup file from a bam file and a reference.
|
|
5
|
|
6 usage: %prog [options]
|
|
7 -i, --input=i: raw snp call file chr:pos
|
|
8 -o, --output1=o: novel snp calls in file
|
|
9 -c, --output2=c: filtered novel SNPs associated with codons
|
|
10 -K, --known_snps=k: known SNPs for filtering (sorted chr:pos file)
|
|
11 -C, --codon=C: codon lookup file (sorted chr:pos)
|
|
12
|
|
13 """
|
|
14
|
|
15 #my $cmd7 = "sort -S 2000M -k 1 $snps | join -a 1 - $known | grep -v dbS | grep -v Vent | grep -v Yor | grep -v Wats | sort -S 2000 -k 1 > $out";
|
|
16 #my $cmd8 = "join $codon $snps\_novel.txt > $snps\_novel." . $base . "codon";
|
|
17
|
|
18 import os, shutil, subprocess, sys, tempfile
|
|
19 from galaxy import eggs
|
|
20 import pkg_resources; pkg_resources.require( "bx-python" )
|
|
21 from bx.cookbook import doc_optparse
|
6
|
22 os.environ['LC_COLLATE'] = 'C'
|
0
|
23
|
|
24 def stop_err( msg ):
|
|
25 sys.stderr.write( '%s\n' % msg )
|
|
26 sys.exit()
|
|
27
|
|
28 def __main__():
|
|
29 #Parse Command Line
|
|
30 options, args = doc_optparse.parse( __doc__ )
|
|
31 # if options.known_snps == "" or options.input == "" or options.codon or "":
|
|
32 # print('Error, required arguments not provided\n')
|
|
33 # return(1)
|
|
34 tmpDir = tempfile.mkdtemp()
|
|
35 #prepare basic filter_snvmix command
|
|
36 filter_cmd = "sort -S 2G -k 1 %s | join -a 1 - %s | grep -v dbS | grep -v Vent | grep -v Yor | grep -v Wats | sort -S 2G -k 1 > %s"
|
|
37 try:
|
|
38 filter_cmd = filter_cmd % ( options.input, options.known_snps, options.output1 )
|
|
39 #run command
|
|
40 #print(filter_cmd)
|
|
41 tmp = tempfile.NamedTemporaryFile( dir=tmpDir ).name
|
|
42 tmp_stderr = open( tmp, 'wb' )
|
|
43 proc = subprocess.Popen( args=filter_cmd, shell=True, cwd=tmpDir, stderr=tmp_stderr.fileno() )
|
|
44 returncode = proc.wait()
|
|
45 tmp_stderr.close()
|
|
46 #did it succeed?
|
|
47 # get stderr, allowing for case where it's very large
|
|
48 tmp_stderr = open( tmp, 'rb' )
|
|
49 stderr = ''
|
|
50 while True:
|
|
51 stderr += tmp_stderr.read( )
|
|
52 if not stderr:
|
|
53 break
|
|
54 tmp_stderr.close()
|
|
55 if returncode != 0:
|
|
56 raise Exception, stderr
|
|
57 except Exception, e:
|
|
58 stop_err( 'Error running filter command\n' + str( e ) )
|
|
59
|
|
60 # check that there are results in the output file
|
|
61 if os.path.getsize( options.output1 ) > 0:
|
|
62 sys.stdout.write( 'wrote output1' )
|
|
63 else:
|
|
64 stop_err( 'The output file is empty. All SNVs might have been known or there may be an error with your input file or settings.' )
|
|
65
|
|
66 codon_cmd = "join %s %s > %s"
|
|
67 try:
|
|
68 codon_cmd = codon_cmd % ( options.codon, options.output1, options.output2 )
|
|
69 #run command
|
|
70 #print(codon_cmd)
|
|
71 tmp = tempfile.NamedTemporaryFile( dir=tmpDir ).name
|
|
72 tmp_stderr = open( tmp, 'wb' )
|
|
73 proc = subprocess.Popen( args=codon_cmd, shell=True, cwd=tmpDir, stderr=tmp_stderr.fileno() )
|
|
74 returncode = proc.wait()
|
|
75 tmp_stderr.close()
|
|
76 #did it succeed?
|
|
77 # get stderr, allowing for case where it's very large
|
|
78 tmp_stderr = open( tmp, 'rb' )
|
|
79 stderr = ''
|
|
80 while True:
|
|
81 stderr += tmp_stderr.read()
|
|
82 if not stderr:
|
|
83 break
|
|
84 tmp_stderr.close()
|
|
85 if returncode != 0:
|
|
86 raise Exception, stderr
|
|
87 except Exception, e:
|
|
88 stop_err( 'Error running codon command\n' + str( e ) )
|
|
89
|
|
90 # check that there are results in the output file
|
|
91 if os.path.getsize( options.output1 ) > 0:
|
|
92 sys.stdout.write( 'wrote output2' )
|
|
93 else:
|
|
94 stop_err( 'The output file is empty. All SNVs might have been intronic or intergenic or there may be an error with your input file or settings.' )
|
|
95
|
|
96
|
|
97
|
|
98 if __name__ == "__main__" : __main__()
|