annotate SMART/DiffExpAnal/tophat_parallel.py @ 18:94ab73e8a190

Uploaded
author m-zytnicki
date Mon, 29 Apr 2013 03:20:15 -0400
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
18
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
1
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
2 #!/usr/bin/env python
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
3
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
4 import optparse, os, shutil, subprocess, sys, tempfile, fileinput, tarfile,random
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
5
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
6 def stop_err( msg ):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
7 sys.stderr.write( "%s\n" % msg )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
8 sys.exit()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
9
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
10 def toTar(tarFileName, accepted_hits_outputNames):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
11 fileName = os.path.splitext(tarFileName)[0]
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
12 fileNameBaseName = os.path.basename(fileName)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
13 dir = os.path.dirname(tarFileName)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
14 tfile = tarfile.open(tarFileName + ".tmp.tar", "w")
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
15 currentPath = os.getcwd()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
16 os.chdir(dir)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
17 for file in accepted_hits_outputNames:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
18 relativeFileName = os.path.basename(file)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
19 tfile.add(relativeFileName)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
20 os.system("mv %s %s" % (tarFileName + ".tmp.tar", tarFileName))
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
21 tfile.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
22 os.chdir(currentPath)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
23
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
24
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
25 def __main__():
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
26 #Parse Command Line
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
27 parser = optparse.OptionParser()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
28 parser.add_option('-o', '--outputTxtFile', dest='outputTxtFile', help='for Differential expression analysis pipeline, new output option gives a txt output containing the list of mapping results.')
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
29 parser.add_option('-t', '--tar', dest='outputTar', default=None, help='output all accepted hits results in a tar file.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
30 parser.add_option( '-p', '--num-threads', dest='num_threads', help='Use this many threads to align reads. The default is 1.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
31 parser.add_option( '-C', '--color-space', dest='color_space', action='store_true', help='This indicates color-space data' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
32 parser.add_option( '-J', '--junctions-output', dest='junctions_output_file', default='junctions_output.bed', help='Junctions output file; formate is BED.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
33 parser.add_option( '-H', '--hits-output', dest='accepted_hits_output_file', default='hits_output_%s.bam' % random.randrange(0, 10000), help='Accepted hits output file; formate is BAM.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
34 parser.add_option( '', '--own-file', dest='own_file', help='' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
35 parser.add_option( '-D', '--indexes-path', dest='index_path', help='Indexes directory; location of .ebwt and .fa files.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
36 parser.add_option( '-r', '--mate-inner-dist', dest='mate_inner_dist', help='This is the expected (mean) inner distance between mate pairs. \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
37 For, example, for paired end runs with fragments selected at 300bp, \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
38 where each end is 50bp, you should set -r to be 200. There is no default, \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
39 and this parameter is required for paired end runs.')
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
40 parser.add_option( '', '--mate-std-dev', dest='mate_std_dev', help='Standard deviation of distribution on inner distances between male pairs.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
41 parser.add_option( '-a', '--min-anchor-length', dest='min_anchor_length',
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
42 help='The "anchor length". TopHat will report junctions spanned by reads with at least this many bases on each side of the junction.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
43 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.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
44 parser.add_option( '-i', '--min-intron-length', dest='min_intron_length',
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
45 help='The minimum intron length. TopHat will ignore donor/acceptor pairs closer than this many bases apart.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
46 parser.add_option( '-I', '--max-intron-length', dest='max_intron_length',
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
47 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.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
48 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)' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
49 parser.add_option( '-g', '--max_multihits', dest='max_multihits', help='Maximum number of alignments to be allowed' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
50 parser.add_option( '', '--initial-read-mismatches', dest='initial_read_mismatches', help='Number of mismatches allowed in the initial read mapping' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
51 parser.add_option( '', '--seg-mismatches', dest='seg_mismatches', help='Number of mismatches allowed in each segment alignment for reads mapped independently' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
52 parser.add_option( '', '--seg-length', dest='seg_length', help='Minimum length of read segments' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
53 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.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
54 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)' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
55 parser.add_option( '', '--max-insertion-length', dest='max_insertion_length', help='The maximum insertion length. The default is 3.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
56 parser.add_option( '', '--max-deletion-length', dest='max_deletion_length', help='The maximum deletion length. The default is 3.' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
57
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
58 # Options for supplying own junctions
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
59 parser.add_option( '-G', '--GTF', dest='gene_model_annotations', help='Supply TopHat with a list of gene model annotations. \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
60 TopHat will use the exon records in this file to build \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
61 a set of known splice junctions for each gene, and will \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
62 attempt to align reads to these junctions even if they \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
63 would not normally be covered by the initial mapping.')
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
64 parser.add_option( '-j', '--raw-juncs', dest='raw_juncs', help='Supply TopHat with a list of raw junctions. Junctions are \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
65 specified one per line, in a tab-delimited format. Records \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
66 look like: <chrom> <left> <right> <+/-> left and right are \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
67 zero-based coordinates, and specify the last character of the \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
68 left sequenced to be spliced to the first character of the right \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
69 sequence, inclusive.')
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
70 parser.add_option( '', '--no-novel-juncs', action="store_true", dest='no_novel_juncs', help="Only look for junctions indicated in the \
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
71 supplied GFF file. (ignored without -G)")
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
72 parser.add_option( '', '--no-novel-indels', action="store_true", dest='no_novel_indels', help="Skip indel search. Indel search is enabled by default.")
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
73 # Types of search.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
74 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.')
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
75 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)')
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
76 parser.add_option( '', '--no-closure-search', action="store_false", dest='closure_search' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
77 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.')
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
78 parser.add_option( '', '--no-coverage-search', action="store_false", dest='coverage_search' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
79 parser.add_option( '', '--min-segment-intron', dest='min_segment_intron', help='Minimum intron length that may be found during split-segment search' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
80 parser.add_option( '', '--max-segment-intron', dest='max_segment_intron', help='Maximum intron length that may be found during split-segment search' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
81 parser.add_option( '', '--min-closure-exon', dest='min_closure_exon', help='Minimum length for exonic hops in potential splice graph' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
82 parser.add_option( '', '--min-closure-intron', dest='min_closure_intron', help='Minimum intron length that may be found during closure search' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
83 parser.add_option( '', '--max-closure-intron', dest='max_closure_intron', help='Maximum intron length that may be found during closure search' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
84 parser.add_option( '', '--min-coverage-intron', dest='min_coverage_intron', help='Minimum intron length that may be found during coverage search' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
85 parser.add_option( '', '--max-coverage-intron', dest='max_coverage_intron', help='Maximum intron length that may be found during coverage search' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
86
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
87 # Wrapper options.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
88 parser.add_option( '-1', '--input1', dest='input1', help='A list of the (forward or single-end) reads files of Sanger FASTQ format, txt format' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
89 #parser.add_option( '-1', '--input1', dest='input1', help='The (forward or single-end) reads file in Sanger FASTQ format' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
90 #parser.add_option( '-2', '--input2', dest='input2', help='The reverse reads file in Sanger FASTQ format' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
91 parser.add_option( '-2', '--input2', dest='input2', help='The list of reverse reads file in Sanger FASTQ format' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
92 parser.add_option( '', '--single-paired', dest='single_paired', help='' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
93 parser.add_option( '', '--settings', dest='settings', help='' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
94
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
95 (options, args) = parser.parse_args()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
96
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
97 # output version # of tool
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
98 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
99 tmp_files = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
100 tmp = tempfile.NamedTemporaryFile().name
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
101 tmp_files.append(tmp)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
102 tmp_stdout = open( tmp, 'wb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
103 proc = subprocess.Popen( args='tophat -v', shell=True, stdout=tmp_stdout )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
104 tmp_stdout.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
105 returncode = proc.wait()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
106 stdout = open( tmp_stdout.name, 'rb' ).readline().strip()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
107 if stdout:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
108 sys.stdout.write( '%s\n' % stdout )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
109 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
110 raise Exception
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
111 except:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
112 sys.stdout.write( 'Could not determine Tophat version\n' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
113
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
114 # Color or base space
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
115 space = ''
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
116 if options.color_space:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
117 space = '-C'
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
118
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
119
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
120 #reads = options.input1
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
121 file = open(options.input1,"r")
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
122 lines = file.readlines()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
123 inputFileNames = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
124 accepted_hits_outputNames = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
125 outputName = options.outputTxtFile
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
126 resDirName = os.path.dirname(outputName) + '/'
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
127 out = open(outputName, "w")
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
128 for line in lines:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
129 tab = line.split()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
130 inputFileNames.append(tab[1])
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
131 aHitOutName = resDirName + tab[0] + '_' + options.accepted_hits_output_file
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
132 accepted_hits_outputNames.append(aHitOutName)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
133 out.write(tab[0] + '\t' + aHitOutName + '\n')
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
134 file.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
135 out.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
136
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
137 if options.input2:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
138 revFile = open(options.input2,"r")
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
139 lines = revFile.readlines()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
140 inputRevFileNames = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
141 for line in lines:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
142 revTab = line.split()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
143 inputRevFileNames.append(revTab[1])
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
144 revFile.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
145
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
146
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
147 # Creat bowtie index if necessary.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
148 tmp_index_dirs = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
149 index_paths = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
150 tmp_index_dir = tempfile.mkdtemp()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
151 tmp_index_dirs.append(tmp_index_dir)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
152 if options.own_file:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
153 index_path = os.path.join( tmp_index_dir, '.'.join( os.path.split( options.own_file )[1].split( '.' )[:-1] ) )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
154 index_paths.append(index_path)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
155 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
156 os.link( options.own_file, index_path + '.fa' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
157 except:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
158 # Tophat prefers (but doesn't require) fasta file to be in same directory, with .fa extension
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
159 pass
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
160 cmd_index = 'bowtie-build %s -f %s %s' % ( space, options.own_file, index_path )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
161 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
162 tmp = tempfile.NamedTemporaryFile( dir=tmp_index_dir ).name
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
163 tmp_stderr = open( tmp, 'wb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
164 proc = subprocess.Popen( args=cmd_index, shell=True, cwd=tmp_index_dir, stderr=tmp_stderr.fileno() )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
165 returncode = proc.wait()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
166 tmp_stderr.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
167 # get stderr, allowing for case where it's very large
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
168 tmp_stderr = open( tmp, 'rb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
169 stderr = ''
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
170 buffsize = 1048576
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
171 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
172 while True:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
173 stderr += tmp_stderr.read( buffsize )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
174 if not stderr or len( stderr ) % buffsize != 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
175 break
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
176 except OverflowError:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
177 pass
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
178 tmp_stderr.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
179 if returncode != 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
180 raise Exception, stderr
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
181 except Exception, e:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
182 if os.path.exists( tmp_index_dir ):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
183 shutil.rmtree( tmp_index_dir )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
184 stop_err( 'Error indexing reference sequence\n' + str( e ) )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
185 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
186 for file in inputFileNames:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
187 tmp_index_dir = tempfile.mkdtemp()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
188 index_path = tmp_index_dir + '/' + os.path.basename(file).split('.')[0]
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
189 index_paths.append(index_path)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
190 tmp_index_dirs.append(tmp_index_dir)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
191
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
192
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
193
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
194 # Build tophat command.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
195 cmds = []
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
196 # for inputFileName in inputFileNames:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
197 for i in range(len(inputFileNames)):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
198 cmd = 'tophat %s %s %s '
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
199 input_files = inputFileNames[i]
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
200 if options.input2:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
201 input_files += ' ' + inputRevFileNames[i]
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
202 opts = '-p %s %s' % ( options.num_threads, space )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
203 if options.single_paired == 'paired':
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
204 opts += '-r %s ' % options.mate_inner_dist
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
205 if options.settings == 'preSet':
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
206 if options.own_file:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
207 cmd = cmd % ( opts, index_paths[0], input_files ) #here add paired end file
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
208 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
209 cmd = cmd % ( opts, index_paths[i], input_files ) #here add paired end file
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
210 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
211 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
212 if int( options.min_anchor_length ) >= 3:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
213 opts += '-a %s ' % options.min_anchor_length
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
214 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
215 raise Exception, 'Minimum anchor length must be 3 or greater'
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
216 opts += '-m %s ' % options.splice_mismatches
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
217 opts += '-i %s ' % options.min_intron_length
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
218 opts += '-I %s ' % options.max_intron_length
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
219 if float( options.junction_filter ) != 0.0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
220 opts += '-F %s ' % options.junction_filter
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
221 opts += '-g %s ' % options.max_multihits
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
222 # Custom junctions options.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
223 if options.gene_model_annotations:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
224 opts += '-G %s ' % options.gene_model_annotations
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
225 if options.raw_juncs:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
226 opts += '-j %s ' % options.raw_juncs
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
227 if options.no_novel_juncs:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
228 opts += '--no-novel-juncs '
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
229 if options.library_type:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
230 opts += '--library-type %s ' % options.library_type
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
231 if options.no_novel_indels:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
232 opts += '--no-novel-indels '
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
233 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
234 if options.max_insertion_length:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
235 opts += '--max-insertion-length %i ' % int( options.max_insertion_length )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
236 if options.max_deletion_length:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
237 opts += '--max-deletion-length %i ' % int( options.max_deletion_length )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
238 # Max options do not work for Tophat v1.2.0, despite documentation to the contrary. (Fixed in version 1.3.1)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
239 # need to warn user of this fact
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
240 #sys.stdout.write( "Max insertion length and max deletion length options don't work in Tophat v1.2.0\n" )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
241
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
242 # Search type options.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
243 if options.coverage_search:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
244 opts += '--coverage-search --min-coverage-intron %s --max-coverage-intron %s ' % ( options.min_coverage_intron, options.max_coverage_intron )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
245 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
246 opts += '--no-coverage-search '
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
247 if options.closure_search:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
248 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 )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
249 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
250 opts += '--no-closure-search '
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
251 if options.microexon_search:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
252 opts += '--microexon-search '
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
253 if options.single_paired == 'paired':
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
254 opts += '--mate-std-dev %s ' % options.mate_std_dev
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
255 if options.initial_read_mismatches:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
256 opts += '--initial-read-mismatches %d ' % int( options.initial_read_mismatches )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
257 if options.seg_mismatches:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
258 opts += '--segment-mismatches %d ' % int( options.seg_mismatches )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
259 if options.seg_length:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
260 opts += '--segment-length %d ' % int( options.seg_length )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
261 if options.min_segment_intron:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
262 opts += '--min-segment-intron %d ' % int( options.min_segment_intron )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
263 if options.max_segment_intron:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
264 opts += '--max-segment-intron %d ' % int( options.max_segment_intron )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
265 if options.own_file:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
266 cmd = cmd % ( opts, index_paths[0], input_files ) #here to add paired end file
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
267 else:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
268 cmd = cmd % ( opts, index_paths[i], input_files ) #here to add paired end file
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
269 except Exception, e:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
270 # Clean up temp dirs
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
271 if os.path.exists( tmp_index_dir ):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
272 shutil.rmtree( tmp_index_dir )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
273 stop_err( 'Something is wrong with the alignment parameters and the alignment could not be run\n' + str( e ) )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
274
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
275 cmds.append(cmd)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
276
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
277 # Run the command line for each file.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
278 for i in range(len(cmds)):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
279 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
280 tmp_out = tempfile.NamedTemporaryFile().name
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
281 tmp_files.append(tmp_out)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
282 tmp_stdout = open( tmp_out, 'wb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
283 tmp_err = tempfile.NamedTemporaryFile().name
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
284 tmp_files.append(tmp_err)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
285 tmp_stderr = open( tmp_err, 'wb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
286 proc = subprocess.Popen( args=cmds[i], shell=True, cwd=".", stdout=tmp_stdout, stderr=tmp_stderr )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
287 returncode = proc.wait()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
288 tmp_stderr.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
289 # get stderr, allowing for case where it's very large
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
290 tmp_stderr = open( tmp_err, 'rb' )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
291 stderr = ''
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
292 buffsize = 1048576
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
293 try:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
294 while True:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
295 stderr += tmp_stderr.read( buffsize )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
296 if not stderr or len( stderr ) % buffsize != 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
297 break
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
298 except OverflowError:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
299 pass
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
300 tmp_stdout.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
301 tmp_stderr.close()
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
302 if returncode != 0:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
303 raise Exception, stderr
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
304
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
305 # Copy output files from tmp directory to specified files.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
306 #shutil.copyfile( os.path.join( "tophat_out", "junctions.bed" ), junctions_outputNames[i] )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
307 shutil.copyfile( os.path.join( "tophat_out", "accepted_hits.bam" ), accepted_hits_outputNames[i] )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
308 # TODO: look for errors in program output.
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
309 except Exception, e:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
310 stop_err( 'Error in tophat:\n' + str( e ) )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
311
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
312 if options.outputTar != None:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
313 toTar(options.outputTar, accepted_hits_outputNames)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
314
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
315
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
316 # Clean up temp dirs
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
317 for tmp_index_dir in tmp_index_dirs:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
318 if os.path.exists( tmp_index_dir ):
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
319 shutil.rmtree( tmp_index_dir )
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
320
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
321 for tmp in tmp_files:
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
322 os.remove(tmp)
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
323
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
324
94ab73e8a190 Uploaded
m-zytnicki
parents:
diff changeset
325 if __name__=="__main__": __main__()