Mercurial > repos > devteam > data_manager_fetch_genome_all_fasta
comparison data_manager/data_manager_fetch_genome_all_fasta.py @ 0:2ebc856bce29 draft
Uploaded data manager definition.
author | devteam |
---|---|
date | Fri, 28 Mar 2014 14:14:24 -0400 |
parents | |
children | ac850912d386 |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:2ebc856bce29 |
---|---|
1 #!/usr/bin/env python | |
2 #Dan Blankenberg | |
3 | |
4 import sys | |
5 import os | |
6 import tempfile | |
7 import shutil | |
8 import optparse | |
9 import urllib2 | |
10 #import uuid | |
11 from ftplib import FTP | |
12 import tarfile | |
13 import zipfile | |
14 import gzip | |
15 import bz2 | |
16 | |
17 from galaxy.util.json import from_json_string, to_json_string | |
18 | |
19 | |
20 CHUNK_SIZE = 2**20 #1mb | |
21 | |
22 def cleanup_before_exit( tmp_dir ): | |
23 if tmp_dir and os.path.exists( tmp_dir ): | |
24 shutil.rmtree( tmp_dir ) | |
25 | |
26 def stop_err(msg): | |
27 sys.stderr.write(msg) | |
28 sys.exit(1) | |
29 | |
30 def get_dbkey_id_name( params, dbkey_description=None ): | |
31 dbkey = params['param_dict']['dbkey'] | |
32 #TODO: ensure sequence_id is unique and does not already appear in location file | |
33 sequence_id = params['param_dict']['sequence_id'] | |
34 if not sequence_id: | |
35 sequence_id = dbkey #uuid.uuid4() generate and use an uuid instead? | |
36 | |
37 sequence_name = params['param_dict']['sequence_name'] | |
38 if not sequence_name: | |
39 sequence_name = dbkey_description | |
40 if not sequence_name: | |
41 sequence_name = dbkey | |
42 return dbkey, sequence_id, sequence_name | |
43 | |
44 def _get_files_in_ftp_path( ftp, path ): | |
45 path_contents = [] | |
46 ftp.retrlines( 'MLSD %s' % ( path ), path_contents.append ) | |
47 return [ line.split( ';' )[ -1 ].lstrip() for line in path_contents ] | |
48 | |
49 def _get_stream_readers_for_tar( file_obj, tmp_dir ): | |
50 fasta_tar = tarfile.open( fileobj=file_obj, mode='r:*' ) | |
51 return [ fasta_tar.extractfile( member ) for member in fasta_tar.getmembers() ] | |
52 | |
53 def _get_stream_readers_for_zip( file_obj, tmp_dir ): | |
54 fasta_zip = zipfile.ZipFile( file_obj, 'r' ) | |
55 rval = [] | |
56 for member in fasta_zip.namelist(): | |
57 fasta_zip.extract( member, tmp_dir ) | |
58 rval.append( open( os.path.join( tmp_dir, member ), 'rb' ) ) | |
59 return rval | |
60 | |
61 def _get_stream_readers_for_gzip( file_obj, tmp_dir ): | |
62 return [ gzip.GzipFile( fileobj=file_obj, mode='rb' ) ] | |
63 | |
64 def _get_stream_readers_for_bz2( file_obj, tmp_dir ): | |
65 return [ bz2.BZ2File( file_obj.name, 'rb' ) ] | |
66 | |
67 def sort_fasta( fasta_filename, sort_method, params ): | |
68 if sort_method is None: | |
69 return | |
70 assert sort_method in SORTING_METHODS, ValueError( "%s is not a valid sorting option." % sort_method ) | |
71 return SORTING_METHODS[ sort_method ]( fasta_filename, params ) | |
72 | |
73 def _move_and_index_fasta_for_sorting( fasta_filename ): | |
74 unsorted_filename = tempfile.NamedTemporaryFile().name | |
75 shutil.move( fasta_filename, unsorted_filename ) | |
76 fasta_offsets = {} | |
77 unsorted_fh = open( unsorted_filename ) | |
78 while True: | |
79 offset = unsorted_fh.tell() | |
80 line = unsorted_fh.readline() | |
81 if not line: | |
82 break | |
83 if line.startswith( ">" ): | |
84 line = line.split( None, 1 )[0][1:] | |
85 fasta_offsets[ line ] = offset | |
86 unsorted_fh.close() | |
87 current_order = map( lambda x: x[1], sorted( map( lambda x: ( x[1], x[0] ), fasta_offsets.items() ) ) ) | |
88 return ( unsorted_filename, fasta_offsets, current_order ) | |
89 | |
90 def _write_sorted_fasta( sorted_names, fasta_offsets, sorted_fasta_filename, unsorted_fasta_filename ): | |
91 unsorted_fh = open( unsorted_fasta_filename ) | |
92 sorted_fh = open( sorted_fasta_filename, 'wb+' ) | |
93 | |
94 for name in sorted_names: | |
95 offset = fasta_offsets[ name ] | |
96 unsorted_fh.seek( offset ) | |
97 sorted_fh.write( unsorted_fh.readline() ) | |
98 while True: | |
99 line = unsorted_fh.readline() | |
100 if not line or line.startswith( ">" ): | |
101 break | |
102 sorted_fh.write( line ) | |
103 unsorted_fh.close() | |
104 sorted_fh.close() | |
105 | |
106 def _sort_fasta_as_is( fasta_filename, params ): | |
107 return | |
108 | |
109 def _sort_fasta_lexicographical( fasta_filename, params ): | |
110 ( unsorted_filename, fasta_offsets, current_order ) = _move_and_index_fasta_for_sorting( fasta_filename ) | |
111 sorted_names = sorted( fasta_offsets.keys() ) | |
112 if sorted_names == current_order: | |
113 shutil.move( unsorted_filename, fasta_filename ) | |
114 else: | |
115 _write_sorted_fasta( sorted_names, fasta_offsets, fasta_filename, unsorted_filename ) | |
116 | |
117 def _sort_fasta_gatk( fasta_filename, params ): | |
118 #This method was added by reviewer request. | |
119 ( unsorted_filename, fasta_offsets, current_order ) = _move_and_index_fasta_for_sorting( fasta_filename ) | |
120 sorted_names = map( str, range( 1, 23 ) ) + [ 'X', 'Y' ] | |
121 #detect if we have chrN, or just N | |
122 has_chr = False | |
123 for chrom in sorted_names: | |
124 if "chr%s" % chrom in current_order: | |
125 has_chr = True | |
126 break | |
127 | |
128 if has_chr: | |
129 sorted_names = map( lambda x: "chr%s" % x, sorted_names) | |
130 sorted_names.insert( 0, "chrM" ) | |
131 else: | |
132 sorted_names.insert( 0, "MT" ) | |
133 sorted_names.extend( map( lambda x: "%s_random" % x, sorted_names ) ) | |
134 | |
135 existing_sorted_names = [] | |
136 for name in sorted_names: | |
137 if name in current_order: | |
138 existing_sorted_names.append( name ) | |
139 for name in current_order: | |
140 #TODO: confirm that non-canonical names do not need to be sorted specially | |
141 if name not in existing_sorted_names: | |
142 existing_sorted_names.append( name ) | |
143 | |
144 if existing_sorted_names == current_order: | |
145 shutil.move( unsorted_filename, fasta_filename ) | |
146 else: | |
147 _write_sorted_fasta( existing_sorted_names, fasta_offsets, fasta_filename, unsorted_filename ) | |
148 | |
149 def _sort_fasta_custom( fasta_filename, params ): | |
150 ( unsorted_filename, fasta_offsets, current_order ) = _move_and_index_fasta_for_sorting( fasta_filename ) | |
151 sorted_names = [] | |
152 for id_repeat in params['param_dict']['sorting']['sequence_identifiers']: | |
153 sorted_names.append( id_repeat[ 'identifier' ] ) | |
154 handle_not_listed = params['param_dict']['sorting']['handle_not_listed']['handle_not_listed_selector'] | |
155 if handle_not_listed.startswith( 'keep' ): | |
156 add_list = [] | |
157 for name in current_order: | |
158 if name not in sorted_names: | |
159 add_list.append( name ) | |
160 if add_list: | |
161 if handle_not_listed == 'keep_append': | |
162 sorted_names.extend( add_list ) | |
163 else: | |
164 add_list.extend( sorted_names ) | |
165 sorted_names = add_list | |
166 if sorted_names == current_order: | |
167 shutil.move( unsorted_filename, fasta_filename ) | |
168 else: | |
169 _write_sorted_fasta( sorted_names, fasta_offsets, fasta_filename, unsorted_filename ) | |
170 | |
171 def download_from_ucsc( data_manager_dict, params, target_directory, dbkey, sequence_id, sequence_name ): | |
172 UCSC_FTP_SERVER = 'hgdownload.cse.ucsc.edu' | |
173 UCSC_CHROM_FA_FILENAME = 'chromFa' | |
174 UCSC_DOWNLOAD_PATH = '/goldenPath/%s/bigZips/' | |
175 COMPRESSED_EXTENSIONS = [ ( '.tar.gz', _get_stream_readers_for_tar ), ( '.tar.bz2', _get_stream_readers_for_tar ), ( '.zip', _get_stream_readers_for_zip ), ( '.fa.gz', _get_stream_readers_for_gzip ), ( '.fa.bz2', _get_stream_readers_for_bz2 ) ] | |
176 | |
177 email = params['param_dict']['__user_email__'] | |
178 if not email: | |
179 email = 'anonymous@example.com' | |
180 | |
181 ucsc_dbkey = params['param_dict']['reference_source']['requested_dbkey'] or dbkey | |
182 ftp = FTP( UCSC_FTP_SERVER ) | |
183 ftp.login( 'anonymous', email ) | |
184 | |
185 ucsc_path = UCSC_DOWNLOAD_PATH % ucsc_dbkey | |
186 path_contents = _get_files_in_ftp_path( ftp, ucsc_path ) | |
187 | |
188 ucsc_file_name = None | |
189 get_stream_reader = None | |
190 ext = None | |
191 for ext, get_stream_reader in COMPRESSED_EXTENSIONS: | |
192 if "%s%s" % ( UCSC_CHROM_FA_FILENAME, ext ) in path_contents: | |
193 ucsc_file_name = "%s%s%s" % ( ucsc_path, UCSC_CHROM_FA_FILENAME, ext ) | |
194 break | |
195 | |
196 if not ucsc_file_name: | |
197 raise Exception( 'Unable to determine filename for UCSC Genome for %s: %s' % ( ucsc_dbkey, path_contents ) ) | |
198 | |
199 | |
200 tmp_dir = tempfile.mkdtemp( prefix='tmp-data-manager-ucsc-' ) | |
201 ucsc_fasta_filename = os.path.join( tmp_dir, "%s%s" % ( UCSC_CHROM_FA_FILENAME, ext ) ) | |
202 | |
203 fasta_base_filename = "%s.fa" % sequence_id | |
204 fasta_filename = os.path.join( target_directory, fasta_base_filename ) | |
205 fasta_writer = open( fasta_filename, 'wb+' ) | |
206 | |
207 tmp_extract_dir = os.path.join ( tmp_dir, 'extracted_fasta' ) | |
208 os.mkdir( tmp_extract_dir ) | |
209 | |
210 tmp_fasta = open( ucsc_fasta_filename, 'wb+' ) | |
211 | |
212 ftp.retrbinary( 'RETR %s' % ucsc_file_name, tmp_fasta.write ) | |
213 | |
214 tmp_fasta.flush() | |
215 tmp_fasta.seek( 0 ) | |
216 | |
217 fasta_readers = get_stream_reader( tmp_fasta, tmp_extract_dir ) | |
218 | |
219 data_table_entry = _stream_fasta_to_file( fasta_readers, target_directory, dbkey, sequence_id, sequence_name, params ) | |
220 _add_data_table_entry( data_manager_dict, data_table_entry ) | |
221 | |
222 for fasta_reader in fasta_readers: | |
223 fasta_reader.close() | |
224 tmp_fasta.close() | |
225 cleanup_before_exit( tmp_dir ) | |
226 | |
227 def download_from_ncbi( data_manager_dict, params, target_directory, dbkey, sequence_id, sequence_name ): | |
228 NCBI_DOWNLOAD_URL = 'http://togows.dbcls.jp/entry/ncbi-nucleotide/%s.fasta' #FIXME: taken from dave's genome manager...why some japan site? | |
229 | |
230 requested_identifier = params['param_dict']['reference_source']['requested_identifier'] | |
231 url = NCBI_DOWNLOAD_URL % requested_identifier | |
232 fasta_reader = urllib2.urlopen( url ) | |
233 | |
234 data_table_entry = _stream_fasta_to_file( fasta_reader, target_directory, dbkey, sequence_id, sequence_name, params ) | |
235 _add_data_table_entry( data_manager_dict, data_table_entry ) | |
236 | |
237 def download_from_url( data_manager_dict, params, target_directory, dbkey, sequence_id, sequence_name ): | |
238 #TODO: we should automatically do decompression here | |
239 urls = filter( bool, map( lambda x: x.strip(), params['param_dict']['reference_source']['user_url'].split( '\n' ) ) ) | |
240 fasta_reader = [ urllib2.urlopen( url ) for url in urls ] | |
241 | |
242 data_table_entry = _stream_fasta_to_file( fasta_reader, target_directory, dbkey, sequence_id, sequence_name, params ) | |
243 _add_data_table_entry( data_manager_dict, data_table_entry ) | |
244 | |
245 def download_from_history( data_manager_dict, params, target_directory, dbkey, sequence_id, sequence_name ): | |
246 #TODO: allow multiple FASTA input files | |
247 input_filename = params['param_dict']['reference_source']['input_fasta'] | |
248 if isinstance( input_filename, list ): | |
249 fasta_reader = [ open( filename, 'rb' ) for filename in input_filename ] | |
250 else: | |
251 fasta_reader = open( input_filename ) | |
252 | |
253 data_table_entry = _stream_fasta_to_file( fasta_reader, target_directory, dbkey, sequence_id, sequence_name, params ) | |
254 _add_data_table_entry( data_manager_dict, data_table_entry ) | |
255 | |
256 def copy_from_directory( data_manager_dict, params, target_directory, dbkey, sequence_id, sequence_name ): | |
257 input_filename = params['param_dict']['reference_source']['fasta_filename'] | |
258 create_symlink = params['param_dict']['reference_source']['create_symlink'] == 'create_symlink' | |
259 if create_symlink: | |
260 data_table_entry = _create_symlink( input_filename, target_directory, dbkey, sequence_id, sequence_name ) | |
261 else: | |
262 if isinstance( input_filename, list ): | |
263 fasta_reader = [ open( filename, 'rb' ) for filename in input_filename ] | |
264 else: | |
265 fasta_reader = open( input_filename ) | |
266 data_table_entry = _stream_fasta_to_file( fasta_reader, target_directory, dbkey, sequence_id, sequence_name, params ) | |
267 _add_data_table_entry( data_manager_dict, data_table_entry ) | |
268 | |
269 def _add_data_table_entry( data_manager_dict, data_table_entry ): | |
270 data_manager_dict['data_tables'] = data_manager_dict.get( 'data_tables', {} ) | |
271 data_manager_dict['data_tables']['all_fasta'] = data_manager_dict['data_tables'].get( 'all_fasta', [] ) | |
272 data_manager_dict['data_tables']['all_fasta'].append( data_table_entry ) | |
273 return data_manager_dict | |
274 | |
275 def _stream_fasta_to_file( fasta_stream, target_directory, dbkey, sequence_id, sequence_name, params, close_stream=True ): | |
276 fasta_base_filename = "%s.fa" % sequence_id | |
277 fasta_filename = os.path.join( target_directory, fasta_base_filename ) | |
278 fasta_writer = open( fasta_filename, 'wb+' ) | |
279 | |
280 if isinstance( fasta_stream, list ) and len( fasta_stream ) == 1: | |
281 fasta_stream = fasta_stream[0] | |
282 | |
283 if isinstance( fasta_stream, list ): | |
284 last_char = None | |
285 for fh in fasta_stream: | |
286 if last_char not in [ None, '\n', '\r' ]: | |
287 fasta_writer.write( '\n' ) | |
288 while True: | |
289 data = fh.read( CHUNK_SIZE ) | |
290 if data: | |
291 fasta_writer.write( data ) | |
292 last_char = data[-1] | |
293 else: | |
294 break | |
295 if close_stream: | |
296 fh.close() | |
297 else: | |
298 while True: | |
299 data = fasta_stream.read( CHUNK_SIZE ) | |
300 if data: | |
301 fasta_writer.write( data ) | |
302 else: | |
303 break | |
304 if close_stream: | |
305 fasta_stream.close() | |
306 | |
307 fasta_writer.close() | |
308 | |
309 sort_fasta( fasta_filename, params['param_dict']['sorting']['sort_selector'], params ) | |
310 | |
311 return dict( value=sequence_id, dbkey=dbkey, name=sequence_name, path=fasta_base_filename ) | |
312 | |
313 def _create_symlink( input_filename, target_directory, dbkey, sequence_id, sequence_name ): | |
314 fasta_base_filename = "%s.fa" % sequence_id | |
315 fasta_filename = os.path.join( target_directory, fasta_base_filename ) | |
316 os.symlink( input_filename, fasta_filename ) | |
317 return dict( value=sequence_id, dbkey=dbkey, name=sequence_name, path=fasta_base_filename ) | |
318 | |
319 | |
320 | |
321 | |
322 REFERENCE_SOURCE_TO_DOWNLOAD = dict( ucsc=download_from_ucsc, ncbi=download_from_ncbi, url=download_from_url, history=download_from_history, directory=copy_from_directory ) | |
323 | |
324 SORTING_METHODS = dict( as_is=_sort_fasta_as_is, lexicographical=_sort_fasta_lexicographical, gatk=_sort_fasta_gatk, custom=_sort_fasta_custom ) | |
325 | |
326 def main(): | |
327 #Parse Command Line | |
328 parser = optparse.OptionParser() | |
329 parser.add_option( '-d', '--dbkey_description', dest='dbkey_description', action='store', type="string", default=None, help='dbkey_description' ) | |
330 (options, args) = parser.parse_args() | |
331 | |
332 filename = args[0] | |
333 | |
334 params = from_json_string( open( filename ).read() ) | |
335 target_directory = params[ 'output_data' ][0]['extra_files_path'] | |
336 os.mkdir( target_directory ) | |
337 data_manager_dict = {} | |
338 | |
339 dbkey, sequence_id, sequence_name = get_dbkey_id_name( params, dbkey_description=options.dbkey_description ) | |
340 | |
341 if dbkey in [ None, '', '?' ]: | |
342 raise Exception( '"%s" is not a valid dbkey. You must specify a valid dbkey.' % ( dbkey ) ) | |
343 | |
344 #Fetch the FASTA | |
345 REFERENCE_SOURCE_TO_DOWNLOAD[ params['param_dict']['reference_source']['reference_source_selector'] ]( data_manager_dict, params, target_directory, dbkey, sequence_id, sequence_name ) | |
346 | |
347 #save info to json file | |
348 open( filename, 'wb' ).write( to_json_string( data_manager_dict ) ) | |
349 | |
350 if __name__ == "__main__": main() |