/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

  1. #!/usr/bin/env python
  2. import os, sys, tempfile
  3. assert sys.version_info[:2] >= (2.4)
  4. def stop_err( msg ):
  5. sys.stderr.write( "%s\n" % msg )
  6. sys.exit()
  7. def check_nib_file( dbkey, GALAXY_DATA_INDEX_DIR ):
  8. nib_file = "%s/alignseq.loc" % GALAXY_DATA_INDEX_DIR
  9. nib_path = ''
  10. nibs = {}
  11. for i, line in enumerate( file( nib_file ) ):
  12. line = line.rstrip( '\r\n' )
  13. if line and not line.startswith( "#" ):
  14. fields = line.split( '\t' )
  15. if len( fields ) < 3:
  16. continue
  17. if fields[0] == 'seq':
  18. nibs[( fields[1] )] = fields[2]
  19. if nibs.has_key( dbkey ):
  20. nib_path = nibs[( dbkey )]
  21. return nib_path
  22. def check_twobit_file( dbkey, GALAXY_DATA_INDEX_DIR ):
  23. twobit_file = "%s/twobit.loc" % GALAXY_DATA_INDEX_DIR
  24. twobit_path = ''
  25. twobits = {}
  26. for i, line in enumerate( file( twobit_file ) ):
  27. line = line.rstrip( '\r\n' )
  28. if line and not line.startswith( "#" ):
  29. fields = line.split( '\t' )
  30. if len( fields ) < 2:
  31. continue
  32. twobits[( fields[0] )] = fields[1]
  33. if twobits.has_key( dbkey ):
  34. twobit_path = twobits[( dbkey )]
  35. return twobit_path
  36. def __main__():
  37. # I/O
  38. source_format = sys.argv[1] # 0: dbkey; 1: upload file
  39. target_file = sys.argv[2]
  40. query_file = sys.argv[3]
  41. output_file = sys.argv[4]
  42. min_iden = sys.argv[5]
  43. tile_size = sys.argv[6]
  44. one_off = sys.argv[7]
  45. try:
  46. float(min_iden)
  47. except:
  48. stop_err('Invalid value for minimal identity.')
  49. try:
  50. test = int(tile_size)
  51. assert test >= 6 and test <= 18
  52. except:
  53. stop_err('Invalid value for tile size. DNA word size must be between 6 and 18.')
  54. try:
  55. test = int(one_off)
  56. assert test >= 0 and test <= int(tile_size)
  57. except:
  58. stop_err('Invalid value for mismatch numbers in the word')
  59. GALAXY_DATA_INDEX_DIR = sys.argv[8]
  60. all_files = []
  61. if source_format == '0':
  62. # check target genome
  63. dbkey = target_file
  64. nib_path = check_nib_file( dbkey, GALAXY_DATA_INDEX_DIR )
  65. twobit_path = check_twobit_file( dbkey, GALAXY_DATA_INDEX_DIR )
  66. if not os.path.exists( nib_path ) and not os.path.exists( twobit_path ):
  67. stop_err("No sequences are available for %s, request them by reporting this error." % dbkey)
  68. # check the query file, see whether all of them are legitimate sequence
  69. if nib_path and os.path.isdir( nib_path ):
  70. compress_files = os.listdir(nib_path)
  71. target_path = nib_path
  72. elif twobit_path:
  73. compress_files = [twobit_path]
  74. target_path = ""
  75. else:
  76. stop_err("Requested genome build has no available sequence.")
  77. for file in compress_files:
  78. file = "%s/%s" % ( target_path, file )
  79. file = os.path.normpath(file)
  80. all_files.append(file)
  81. else:
  82. all_files = [target_file]
  83. for detail_file_path in all_files:
  84. output_tempfile = tempfile.NamedTemporaryFile().name
  85. 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 )
  86. os.system( command )
  87. os.system( 'cat %s >> %s' % ( output_tempfile, output_file ) )
  88. os.remove( output_tempfile )
  89. if __name__ == '__main__': __main__()