|
0
|
1 #!/usr/bin/env python
|
|
|
2
|
|
|
3 import optparse, os, shutil, subprocess, sys, tempfile, fileinput
|
|
|
4
|
|
|
5 def stop_err( msg ):
|
|
|
6 sys.stderr.write( "%s\n" % msg )
|
|
|
7 sys.exit()
|
|
|
8
|
|
|
9 def __main__():
|
|
|
10 #Parse Command Line
|
|
|
11 parser = optparse.OptionParser()
|
|
|
12 parser.add_option( '-p', '--num-threads', dest='num_threads', help='Use this many threads to align reads. The default is 1.' )
|
|
|
13 parser.add_option( '-C', '--color-space', dest='color_space', action='store_true', help='This indicates color-space data' )
|
|
|
14 parser.add_option( '-J', '--junctions-output', dest='junctions_output_file', help='Junctions output file; formate is BED.' )
|
|
|
15 parser.add_option( '-H', '--hits-output', dest='accepted_hits_output_file', help='Accepted hits output file; formate is BAM.' )
|
|
|
16 parser.add_option( '', '--own-file', dest='own_file', help='' )
|
|
|
17 parser.add_option( '-D', '--indexes-path', dest='index_path', help='Indexes directory; location of .ebwt and .fa files.' )
|
|
|
18 parser.add_option( '-r', '--mate-inner-dist', dest='mate_inner_dist', help='This is the expected (mean) inner distance between mate pairs. \
|
|
|
19 For, example, for paired end runs with fragments selected at 300bp, \
|
|
|
20 where each end is 50bp, you should set -r to be 200. There is no default, \
|
|
|
21 and this parameter is required for paired end runs.')
|
|
|
22 parser.add_option( '', '--mate-std-dev', dest='mate_std_dev', help='Standard deviation of distribution on inner distances between male pairs.' )
|
|
|
23 parser.add_option( '-a', '--min-anchor-length', dest='min_anchor_length',
|
|
|
24 help='The "anchor length". TopHat will report junctions spanned by reads with at least this many bases on each side of the junction.' )
|
|
|
25 parser.add_option( '-m', '--splice-mismatches', dest='splice_mismatches', help='The maximum number of mismatches that can appear in the anchor region of a spliced alignment.' )
|
|
|
26 parser.add_option( '-i', '--min-intron-length', dest='min_intron_length',
|
|
|
27 help='The minimum intron length. TopHat will ignore donor/acceptor pairs closer than this many bases apart.' )
|
|
|
28 parser.add_option( '-I', '--max-intron-length', dest='max_intron_length',
|
|
|
29 help='The maximum intron length. When searching for junctions ab initio, TopHat will ignore donor/acceptor pairs farther than this many bases apart, except when such a pair is supported by a split segment alignment of a long read.' )
|
|
|
30 parser.add_option( '-F', '--junction_filter', dest='junction_filter', help='Filter out junctions supported by too few alignments (number of reads divided by average depth of coverage)' )
|
|
|
31 parser.add_option( '-g', '--max_multihits', dest='max_multihits', help='Maximum number of alignments to be allowed' )
|
|
|
32 parser.add_option( '', '--seg-mismatches', dest='seg_mismatches', help='Number of mismatches allowed in each segment alignment for reads mapped independently' )
|
|
|
33 parser.add_option( '', '--seg-length', dest='seg_length', help='Minimum length of read segments' )
|
|
|
34 parser.add_option( '', '--library-type', dest='library_type', help='TopHat will treat the reads as strand specific. Every read alignment will have an XS attribute tag. Consider supplying library type options below to select the correct RNA-seq protocol.' )
|
|
|
35 parser.add_option( '', '--allow-indels', action="store_true", help='Allow indel search. Indel search is disabled by default.' )
|
|
|
36 parser.add_option( '', '--max-insertion-length', dest='max_insertion_length', help='The maximum insertion length. The default is 3.' )
|
|
|
37 parser.add_option( '', '--max-deletion-length', dest='max_deletion_length', help='The maximum deletion length. The default is 3.' )
|
|
|
38
|
|
|
39 # Options for supplying own junctions
|
|
|
40 parser.add_option( '-G', '--GTF', dest='gene_model_annotations', help='Supply TopHat with a list of gene model annotations. \
|
|
|
41 TopHat will use the exon records in this file to build \
|
|
|
42 a set of known splice junctions for each gene, and will \
|
|
|
43 attempt to align reads to these junctions even if they \
|
|
|
44 would not normally be covered by the initial mapping.')
|
|
|
45 parser.add_option( '-j', '--raw-juncs', dest='raw_juncs', help='Supply TopHat with a list of raw junctions. Junctions are \
|
|
|
46 specified one per line, in a tab-delimited format. Records \
|
|
|
47 look like: <chrom> <left> <right> <+/-> left and right are \
|
|
|
48 zero-based coordinates, and specify the last character of the \
|
|
|
49 left sequenced to be spliced to the first character of the right \
|
|
|
50 sequence, inclusive.')
|
|
|
51 parser.add_option( '', '--no-novel-juncs', action="store_true", dest='no_novel_juncs', help="Only look for junctions indicated in the \
|
|
|
52 supplied GFF file. (ignored without -G)")
|
|
|
53 # Types of search.
|
|
|
54 parser.add_option( '', '--microexon-search', action="store_true", dest='microexon_search', help='With this option, the pipeline will attempt to find alignments incident to microexons. Works only for reads 50bp or longer.')
|
|
|
55 parser.add_option( '', '--closure-search', action="store_true", dest='closure_search', help='Enables the mate pair closure-based search for junctions. Closure-based search should only be used when the expected inner distance between mates is small (<= 50bp)')
|
|
|
56 parser.add_option( '', '--no-closure-search', action="store_false", dest='closure_search' )
|
|
|
57 parser.add_option( '', '--coverage-search', action="store_true", dest='coverage_search', help='Enables the coverage based search for junctions. Use when coverage search is disabled by default (such as for reads 75bp or longer), for maximum sensitivity.')
|
|
|
58 parser.add_option( '', '--no-coverage-search', action="store_false", dest='coverage_search' )
|
|
|
59 parser.add_option( '', '--min-segment-intron', dest='min_segment_intron', help='Minimum intron length that may be found during split-segment search' )
|
|
|
60 parser.add_option( '', '--max-segment-intron', dest='max_segment_intron', help='Maximum intron length that may be found during split-segment search' )
|
|
|
61 parser.add_option( '', '--min-closure-exon', dest='min_closure_exon', help='Minimum length for exonic hops in potential splice graph' )
|
|
|
62 parser.add_option( '', '--min-closure-intron', dest='min_closure_intron', help='Minimum intron length that may be found during closure search' )
|
|
|
63 parser.add_option( '', '--max-closure-intron', dest='max_closure_intron', help='Maximum intron length that may be found during closure search' )
|
|
|
64 parser.add_option( '', '--min-coverage-intron', dest='min_coverage_intron', help='Minimum intron length that may be found during coverage search' )
|
|
|
65 parser.add_option( '', '--max-coverage-intron', dest='max_coverage_intron', help='Maximum intron length that may be found during coverage search' )
|
|
|
66
|
|
|
67 # Wrapper options.
|
|
|
68 parser.add_option( '-1', '--input1', dest='input1', help='The (forward or single-end) reads file in Sanger FASTQ format' )
|
|
|
69 parser.add_option( '-2', '--input2', dest='input2', help='The reverse reads file in Sanger FASTQ format' )
|
|
|
70 parser.add_option( '', '--single-paired', dest='single_paired', help='' )
|
|
|
71 parser.add_option( '', '--settings', dest='settings', help='' )
|
|
|
72
|
|
|
73 (options, args) = parser.parse_args()
|
|
|
74
|
|
|
75 # output version # of tool
|
|
|
76 try:
|
|
|
77 tmp = tempfile.NamedTemporaryFile().name
|
|
|
78 tmp_stdout = open( tmp, 'wb' )
|
|
|
79 proc = subprocess.Popen( args='tophat -v', shell=True, stdout=tmp_stdout )
|
|
|
80 tmp_stdout.close()
|
|
|
81 returncode = proc.wait()
|
|
|
82 stdout = open( tmp_stdout.name, 'rb' ).readline().strip()
|
|
|
83 if stdout:
|
|
|
84 sys.stdout.write( '%s\n' % stdout )
|
|
|
85 else:
|
|
|
86 raise Exception
|
|
|
87 except:
|
|
|
88 sys.stdout.write( 'Could not determine Tophat version\n' )
|
|
|
89
|
|
|
90 # Color or base space
|
|
|
91 space = ''
|
|
|
92 if options.color_space:
|
|
|
93 space = '-C'
|
|
|
94
|
|
|
95 # Creat bowtie index if necessary.
|
|
|
96 tmp_index_dir = tempfile.mkdtemp()
|
|
|
97 if options.own_file:
|
|
|
98 index_path = os.path.join( tmp_index_dir, '.'.join( os.path.split( options.own_file )[1].split( '.' )[:-1] ) )
|
|
|
99 try:
|
|
|
100 os.link( options.own_file, index_path + '.fa' )
|
|
|
101 except:
|
|
|
102 # Tophat prefers (but doesn't require) fasta file to be in same directory, with .fa extension
|
|
|
103 pass
|
|
|
104 cmd_index = 'bowtie-build %s -f %s %s' % ( space, options.own_file, index_path )
|
|
|
105 try:
|
|
|
106 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
|
|
|
107 tmp_stderr = open( tmp, 'wb' )
|
|
|
108 proc = subprocess.Popen( args=cmd_index, shell=True, cwd=tmp_index_dir, stderr=tmp_stderr.fileno() )
|
|
|
109 returncode = proc.wait()
|
|
|
110 tmp_stderr.close()
|
|
|
111 # get stderr, allowing for case where it's very large
|
|
|
112 tmp_stderr = open( tmp, 'rb' )
|
|
|
113 stderr = ''
|
|
|
114 buffsize = 1048576
|
|
|
115 try:
|
|
|
116 while True:
|
|
|
117 stderr += tmp_stderr.read( buffsize )
|
|
|
118 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
119 break
|
|
|
120 except OverflowError:
|
|
|
121 pass
|
|
|
122 tmp_stderr.close()
|
|
|
123 if returncode != 0:
|
|
|
124 raise Exception, stderr
|
|
|
125 except Exception, e:
|
|
|
126 if os.path.exists( tmp_index_dir ):
|
|
|
127 shutil.rmtree( tmp_index_dir )
|
|
|
128 stop_err( 'Error indexing reference sequence\n' + str( e ) )
|
|
|
129 else:
|
|
|
130 index_path = options.index_path
|
|
|
131
|
|
|
132 # Build tophat command.
|
|
|
133 cmd = 'tophat %s %s %s'
|
|
|
134 reads = options.input1
|
|
|
135 if options.input2:
|
|
|
136 reads += ' ' + options.input2
|
|
|
137 opts = '-p %s %s' % ( options.num_threads, space )
|
|
|
138 if options.single_paired == 'paired':
|
|
|
139 opts += ' -r %s' % options.mate_inner_dist
|
|
|
140 if options.settings == 'preSet':
|
|
|
141 cmd = cmd % ( opts, index_path, reads )
|
|
|
142 else:
|
|
|
143 try:
|
|
|
144 if int( options.min_anchor_length ) >= 3:
|
|
|
145 opts += ' -a %s' % options.min_anchor_length
|
|
|
146 else:
|
|
|
147 raise Exception, 'Minimum anchor length must be 3 or greater'
|
|
|
148 opts += ' -m %s' % options.splice_mismatches
|
|
|
149 opts += ' -i %s' % options.min_intron_length
|
|
|
150 opts += ' -I %s' % options.max_intron_length
|
|
|
151 if float( options.junction_filter ) != 0.0:
|
|
|
152 opts += ' -F %s' % options.junction_filter
|
|
|
153 opts += ' -g %s' % options.max_multihits
|
|
|
154 # Custom junctions options.
|
|
|
155 if options.gene_model_annotations:
|
|
|
156 opts += ' -G %s' % options.gene_model_annotations
|
|
|
157 if options.raw_juncs:
|
|
|
158 opts += ' -j %s' % options.raw_juncs
|
|
|
159 if options.no_novel_juncs:
|
|
|
160 opts += ' --no-novel-juncs'
|
|
|
161 if options.library_type:
|
|
|
162 opts += ' --library-type %s' % options.library_type
|
|
|
163 if options.allow_indels:
|
|
|
164 # Max options do not work for Tophat v1.2.0, despite documentation to the contrary.
|
|
|
165 opts += ' --allow-indels'
|
|
|
166 #opts += ' --max-insertion-length %i --max-deletion-length %i' % ( int( options.max_insertion_length ), int( options.max_deletion_length ) )
|
|
|
167 # need to warn user of this fact
|
|
|
168 sys.stdout.write( "Max insertion length and max deletion length options don't work in Tophat v1.2.0\n" )
|
|
|
169
|
|
|
170 # Search type options.
|
|
|
171 if options.coverage_search:
|
|
|
172 opts += ' --coverage-search --min-coverage-intron %s --max-coverage-intron %s' % ( options.min_coverage_intron, options.max_coverage_intron )
|
|
|
173 else:
|
|
|
174 opts += ' --no-coverage-search'
|
|
|
175 if options.closure_search:
|
|
|
176 opts += ' --closure-search --min-closure-exon %s --min-closure-intron %s --max-closure-intron %s' % ( options.min_closure_exon, options.min_closure_intron, options.max_closure_intron )
|
|
|
177 else:
|
|
|
178 opts += ' --no-closure-search'
|
|
|
179 if options.microexon_search:
|
|
|
180 opts += ' --microexon-search'
|
|
|
181 if options.single_paired == 'paired':
|
|
|
182 opts += ' --mate-std-dev %s' % options.mate_std_dev
|
|
|
183 if options.seg_mismatches:
|
|
|
184 opts += ' --segment-mismatches %d' % int( options.seg_mismatches )
|
|
|
185 if options.seg_length:
|
|
|
186 opts += ' --segment-length %d' % int( options.seg_length )
|
|
|
187 if options.min_segment_intron:
|
|
|
188 opts += ' --min-segment-intron %d' % int( options.min_segment_intron )
|
|
|
189 if options.max_segment_intron:
|
|
|
190 opts += ' --max-segment-intron %d' % int( options.max_segment_intron )
|
|
|
191 cmd = cmd % ( opts, index_path, reads )
|
|
|
192 except Exception, e:
|
|
|
193 # Clean up temp dirs
|
|
|
194 if os.path.exists( tmp_index_dir ):
|
|
|
195 shutil.rmtree( tmp_index_dir )
|
|
|
196 stop_err( 'Something is wrong with the alignment parameters and the alignment could not be run\n' + str( e ) )
|
|
|
197 #print cmd
|
|
|
198
|
|
|
199 # Run
|
|
|
200 try:
|
|
|
201 tmp_out = tempfile.NamedTemporaryFile().name
|
|
|
202 tmp_stdout = open( tmp_out, 'wb' )
|
|
|
203 tmp_err = tempfile.NamedTemporaryFile().name
|
|
|
204 tmp_stderr = open( tmp_err, 'wb' )
|
|
|
205 proc = subprocess.Popen( args=cmd, shell=True, cwd=".", stdout=tmp_stdout, stderr=tmp_stderr )
|
|
|
206 returncode = proc.wait()
|
|
|
207 tmp_stderr.close()
|
|
|
208 # get stderr, allowing for case where it's very large
|
|
|
209 tmp_stderr = open( tmp_err, 'rb' )
|
|
|
210 stderr = ''
|
|
|
211 buffsize = 1048576
|
|
|
212 try:
|
|
|
213 while True:
|
|
|
214 stderr += tmp_stderr.read( buffsize )
|
|
|
215 if not stderr or len( stderr ) % buffsize != 0:
|
|
|
216 break
|
|
|
217 except OverflowError:
|
|
|
218 pass
|
|
|
219 tmp_stdout.close()
|
|
|
220 tmp_stderr.close()
|
|
|
221 if returncode != 0:
|
|
|
222 raise Exception, stderr
|
|
|
223
|
|
|
224 # Copy output files from tmp directory to specified files.
|
|
|
225 shutil.copyfile( os.path.join( "tophat_out", "junctions.bed" ), options.junctions_output_file )
|
|
|
226 shutil.copyfile( os.path.join( "tophat_out", "accepted_hits.bam" ), options.accepted_hits_output_file )
|
|
|
227
|
|
|
228 # TODO: look for errors in program output.
|
|
|
229 except Exception, e:
|
|
|
230 stop_err( 'Error in tophat:\n' + str( e ) )
|
|
|
231
|
|
|
232 # Clean up temp dirs
|
|
|
233 if os.path.exists( tmp_index_dir ):
|
|
|
234 shutil.rmtree( tmp_index_dir )
|
|
|
235
|
|
|
236 if __name__=="__main__": __main__()
|
