annotate tophat_wrapper.py @ 3:2ad64c5bb5f4 draft default tip

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