Mercurial > repos > ryanmorin > nextgen_variant_identification
comparison SNV/snp_filters.py @ 0:74f5ea818cea
Uploaded
author | ryanmorin |
---|---|
date | Wed, 12 Oct 2011 19:50:38 -0400 |
parents | |
children | 361d6506850a |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:74f5ea818cea |
---|---|
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 | |
22 | |
23 def stop_err( msg ): | |
24 sys.stderr.write( '%s\n' % msg ) | |
25 sys.exit() | |
26 | |
27 def __main__(): | |
28 #Parse Command Line | |
29 options, args = doc_optparse.parse( __doc__ ) | |
30 # if options.known_snps == "" or options.input == "" or options.codon or "": | |
31 # print('Error, required arguments not provided\n') | |
32 # return(1) | |
33 tmpDir = tempfile.mkdtemp() | |
34 #prepare basic filter_snvmix command | |
35 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" | |
36 try: | |
37 filter_cmd = filter_cmd % ( options.input, options.known_snps, options.output1 ) | |
38 #run command | |
39 #print(filter_cmd) | |
40 tmp = tempfile.NamedTemporaryFile( dir=tmpDir ).name | |
41 tmp_stderr = open( tmp, 'wb' ) | |
42 proc = subprocess.Popen( args=filter_cmd, shell=True, cwd=tmpDir, stderr=tmp_stderr.fileno() ) | |
43 returncode = proc.wait() | |
44 tmp_stderr.close() | |
45 #did it succeed? | |
46 # get stderr, allowing for case where it's very large | |
47 tmp_stderr = open( tmp, 'rb' ) | |
48 stderr = '' | |
49 while True: | |
50 stderr += tmp_stderr.read( ) | |
51 if not stderr: | |
52 break | |
53 tmp_stderr.close() | |
54 if returncode != 0: | |
55 raise Exception, stderr | |
56 except Exception, e: | |
57 stop_err( 'Error running filter command\n' + str( e ) ) | |
58 | |
59 # check that there are results in the output file | |
60 if os.path.getsize( options.output1 ) > 0: | |
61 sys.stdout.write( 'wrote output1' ) | |
62 else: | |
63 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.' ) | |
64 | |
65 codon_cmd = "join %s %s > %s" | |
66 try: | |
67 codon_cmd = codon_cmd % ( options.codon, options.output1, options.output2 ) | |
68 #run command | |
69 #print(codon_cmd) | |
70 tmp = tempfile.NamedTemporaryFile( dir=tmpDir ).name | |
71 tmp_stderr = open( tmp, 'wb' ) | |
72 proc = subprocess.Popen( args=codon_cmd, shell=True, cwd=tmpDir, stderr=tmp_stderr.fileno() ) | |
73 returncode = proc.wait() | |
74 tmp_stderr.close() | |
75 #did it succeed? | |
76 # get stderr, allowing for case where it's very large | |
77 tmp_stderr = open( tmp, 'rb' ) | |
78 stderr = '' | |
79 while True: | |
80 stderr += tmp_stderr.read() | |
81 if not stderr: | |
82 break | |
83 tmp_stderr.close() | |
84 if returncode != 0: | |
85 raise Exception, stderr | |
86 except Exception, e: | |
87 stop_err( 'Error running codon command\n' + str( e ) ) | |
88 | |
89 # check that there are results in the output file | |
90 if os.path.getsize( options.output1 ) > 0: | |
91 sys.stdout.write( 'wrote output2' ) | |
92 else: | |
93 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.' ) | |
94 | |
95 | |
96 | |
97 if __name__ == "__main__" : __main__() |