annotate sam_pileup.py @ 3:890d97772e2a draft

Uploaded
author devteam
date Thu, 09 Jan 2014 14:28:39 -0500
parents e7d863c5c5d6
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
1 #!/usr/bin/env python
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
2
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
3 """
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
4 Creates a pileup file from a bam file and a reference.
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
5
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
6 usage: %prog [options]
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
7 -p, --input1=p: bam file
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
8 -o, --output1=o: Output pileup
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
9 -R, --ref=R: Reference file type
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
10 -n, --ownFile=n: User-supplied fasta reference file
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
11 -b, --bamIndex=b: BAM index file
1
e7d863c5c5d6 Update sam_pileup to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
12 -g, --index=g: Path of the indexed reference genome
0
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
13 -s, --lastCol=s: Print the mapping quality as the last column
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
14 -i, --indels=i: Only output lines containing indels
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
15 -M, --mapCap=M: Cap mapping quality
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
16 -c, --consensus=c: Call the consensus sequence using MAQ consensu model
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
17 -T, --theta=T: Theta paramter (error dependency coefficient)
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
18 -N, --hapNum=N: Number of haplotypes in sample
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
19 -r, --fraction=r: Expected fraction of differences between a pair of haplotypes
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
20 -I, --phredProb=I: Phred probability of an indel in sequencing/prep
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
21
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
22 """
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
23
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
24 import os, shutil, subprocess, sys, tempfile
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
25 from galaxy import eggs
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
26 import pkg_resources; pkg_resources.require( "bx-python" )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
27 from bx.cookbook import doc_optparse
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
28
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
29 def stop_err( msg ):
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
30 sys.stderr.write( '%s\n' % msg )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
31 sys.exit()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
32
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
33 def __main__():
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
34 #Parse Command Line
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
35 options, args = doc_optparse.parse( __doc__ )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
36 # output version # of tool
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
37 try:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
38 tmp = tempfile.NamedTemporaryFile().name
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
39 tmp_stdout = open( tmp, 'wb' )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
40 proc = subprocess.Popen( args='samtools 2>&1', shell=True, stdout=tmp_stdout )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
41 tmp_stdout.close()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
42 returncode = proc.wait()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
43 stdout = None
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
44 for line in open( tmp_stdout.name, 'rb' ):
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
45 if line.lower().find( 'version' ) >= 0:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
46 stdout = line.strip()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
47 break
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
48 if stdout:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
49 sys.stdout.write( 'Samtools %s\n' % stdout )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
50 else:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
51 raise Exception
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
52 except:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
53 sys.stdout.write( 'Could not determine Samtools version\n' )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
54 #prepare file names
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
55 tmpDir = tempfile.mkdtemp()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
56 tmpf0 = tempfile.NamedTemporaryFile( dir=tmpDir )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
57 tmpf0_name = tmpf0.name
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
58 tmpf0.close()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
59 tmpf0bam_name = '%s.bam' % tmpf0_name
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
60 tmpf0bambai_name = '%s.bam.bai' % tmpf0_name
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
61 tmpf1 = tempfile.NamedTemporaryFile( dir=tmpDir )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
62 tmpf1_name = tmpf1.name
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
63 tmpf1.close()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
64 #link bam and bam index to working directory (can't move because need to leave original)
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
65 os.symlink( options.input1, tmpf0bam_name )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
66 os.symlink( options.bamIndex, tmpf0bambai_name )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
67 #get parameters for pileup command
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
68 if options.lastCol == 'yes':
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
69 lastCol = '-s'
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
70 else:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
71 lastCol = ''
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
72 if options.indels == 'yes':
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
73 indels = '-i'
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
74 else:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
75 indels = ''
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
76 opts = '%s %s -M %s' % ( lastCol, indels, options.mapCap )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
77 if options.consensus == 'yes':
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
78 opts += ' -c -T %s -N %s -r %s -I %s' % ( options.theta, options.hapNum, options.fraction, options.phredProb )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
79 #prepare basic pileup command
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
80 cmd = 'samtools pileup %s -f %s %s > %s'
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
81 try:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
82 # have to nest try-except in try-finally to handle 2.4
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
83 try:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
84 #index reference if necessary and prepare pileup command
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
85 if options.ref == 'indexed':
1
e7d863c5c5d6 Update sam_pileup to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
86 if not os.path.exists( "%s.fai" % options.index ):
e7d863c5c5d6 Update sam_pileup to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
87 raise Exception, "Indexed genome %s not present, request it by reporting this error." % options.index
e7d863c5c5d6 Update sam_pileup to use the fasta_indexes data table.
Dave Bouvier <dave@bx.psu.edu>
parents: 0
diff changeset
88 cmd = cmd % ( opts, options.index, tmpf0bam_name, options.output1 )
0
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
89 elif options.ref == 'history':
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
90 os.symlink( options.ownFile, tmpf1_name )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
91 cmdIndex = 'samtools faidx %s' % ( tmpf1_name )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
92 tmp = tempfile.NamedTemporaryFile( dir=tmpDir ).name
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
93 tmp_stderr = open( tmp, 'wb' )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
94 proc = subprocess.Popen( args=cmdIndex, shell=True, cwd=tmpDir, stderr=tmp_stderr.fileno() )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
95 returncode = proc.wait()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
96 tmp_stderr.close()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
97 # get stderr, allowing for case where it's very large
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
98 tmp_stderr = open( tmp, 'rb' )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
99 stderr = ''
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
100 buffsize = 1048576
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
101 try:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
102 while True:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
103 stderr += tmp_stderr.read( buffsize )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
104 if not stderr or len( stderr ) % buffsize != 0:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
105 break
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
106 except OverflowError:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
107 pass
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
108 tmp_stderr.close()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
109 #did index succeed?
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
110 if returncode != 0:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
111 raise Exception, 'Error creating index file\n' + stderr
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
112 cmd = cmd % ( opts, tmpf1_name, tmpf0bam_name, options.output1 )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
113 #perform pileup command
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
114 tmp = tempfile.NamedTemporaryFile( dir=tmpDir ).name
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
115 tmp_stderr = open( tmp, 'wb' )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
116 proc = subprocess.Popen( args=cmd, shell=True, cwd=tmpDir, stderr=tmp_stderr.fileno() )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
117 returncode = proc.wait()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
118 tmp_stderr.close()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
119 #did it succeed?
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
120 # get stderr, allowing for case where it's very large
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
121 tmp_stderr = open( tmp, 'rb' )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
122 stderr = ''
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
123 buffsize = 1048576
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
124 try:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
125 while True:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
126 stderr += tmp_stderr.read( buffsize )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
127 if not stderr or len( stderr ) % buffsize != 0:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
128 break
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
129 except OverflowError:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
130 pass
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
131 tmp_stderr.close()
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
132 if returncode != 0:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
133 raise Exception, stderr
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
134 except Exception, e:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
135 stop_err( 'Error running Samtools pileup tool\n' + str( e ) )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
136 finally:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
137 #clean up temp files
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
138 if os.path.exists( tmpDir ):
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
139 shutil.rmtree( tmpDir )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
140 # check that there are results in the output file
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
141 if os.path.getsize( options.output1 ) > 0:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
142 sys.stdout.write( 'Converted BAM to pileup' )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
143 else:
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
144 stop_err( 'The output file is empty. Your input file may have had no matches, or there may be an error with your input file or settings.' )
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
145
95612c159681 Uploaded tool tarball.
devteam
parents:
diff changeset
146 if __name__ == "__main__" : __main__()