/tools/metag_tools/blat_wrapper.py
https://bitbucket.org/cistrome/cistrome-harvard/ · Python · 104 lines · 85 code · 15 blank · 4 comment · 27 complexity · f8ce5202f15f392e643d20f4c855bce6 MD5 · raw file
- #!/usr/bin/env python
- import os, sys, tempfile
- assert sys.version_info[:2] >= (2.4)
- def stop_err( msg ):
- sys.stderr.write( "%s\n" % msg )
- sys.exit()
-
- def check_nib_file( dbkey, GALAXY_DATA_INDEX_DIR ):
- nib_file = "%s/alignseq.loc" % GALAXY_DATA_INDEX_DIR
- nib_path = ''
- nibs = {}
- for i, line in enumerate( file( nib_file ) ):
- line = line.rstrip( '\r\n' )
- if line and not line.startswith( "#" ):
- fields = line.split( '\t' )
- if len( fields ) < 3:
- continue
- if fields[0] == 'seq':
- nibs[( fields[1] )] = fields[2]
- if nibs.has_key( dbkey ):
- nib_path = nibs[( dbkey )]
- return nib_path
- def check_twobit_file( dbkey, GALAXY_DATA_INDEX_DIR ):
- twobit_file = "%s/twobit.loc" % GALAXY_DATA_INDEX_DIR
- twobit_path = ''
- twobits = {}
- for i, line in enumerate( file( twobit_file ) ):
- line = line.rstrip( '\r\n' )
- if line and not line.startswith( "#" ):
- fields = line.split( '\t' )
- if len( fields ) < 2:
- continue
- twobits[( fields[0] )] = fields[1]
- if twobits.has_key( dbkey ):
- twobit_path = twobits[( dbkey )]
- return twobit_path
- def __main__():
- # I/O
- source_format = sys.argv[1] # 0: dbkey; 1: upload file
- target_file = sys.argv[2]
- query_file = sys.argv[3]
- output_file = sys.argv[4]
- min_iden = sys.argv[5]
- tile_size = sys.argv[6]
- one_off = sys.argv[7]
-
- try:
- float(min_iden)
- except:
- stop_err('Invalid value for minimal identity.')
-
- try:
- test = int(tile_size)
- assert test >= 6 and test <= 18
- except:
- stop_err('Invalid value for tile size. DNA word size must be between 6 and 18.')
-
- try:
- test = int(one_off)
- assert test >= 0 and test <= int(tile_size)
- except:
- stop_err('Invalid value for mismatch numbers in the word')
-
- GALAXY_DATA_INDEX_DIR = sys.argv[8]
- all_files = []
- if source_format == '0':
- # check target genome
- dbkey = target_file
- nib_path = check_nib_file( dbkey, GALAXY_DATA_INDEX_DIR )
- twobit_path = check_twobit_file( dbkey, GALAXY_DATA_INDEX_DIR )
- if not os.path.exists( nib_path ) and not os.path.exists( twobit_path ):
- stop_err("No sequences are available for %s, request them by reporting this error." % dbkey)
-
- # check the query file, see whether all of them are legitimate sequence
- if nib_path and os.path.isdir( nib_path ):
- compress_files = os.listdir(nib_path)
- target_path = nib_path
- elif twobit_path:
- compress_files = [twobit_path]
- target_path = ""
- else:
- stop_err("Requested genome build has no available sequence.")
-
- for file in compress_files:
- file = "%s/%s" % ( target_path, file )
- file = os.path.normpath(file)
- all_files.append(file)
- else:
- all_files = [target_file]
-
- for detail_file_path in all_files:
- output_tempfile = tempfile.NamedTemporaryFile().name
- command = "blat %s %s %s -oneOff=%s -tileSize=%s -minIdentity=%s -mask=lower -noHead -out=pslx 2>&1" % ( detail_file_path, query_file, output_tempfile, one_off, tile_size, min_iden )
- os.system( command )
- os.system( 'cat %s >> %s' % ( output_tempfile, output_file ) )
- os.remove( output_tempfile )
-
- if __name__ == '__main__': __main__()