/tools/metag_tools/blat_wrapper.py

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