[2] | 1 | #!/usr/bin/env python |
---|
| 2 | #Guruprasad Ananda |
---|
| 3 | """ |
---|
| 4 | Converts coordinates from one build/assembly to another using liftOver binary and mapping files downloaded from UCSC. |
---|
| 5 | """ |
---|
| 6 | |
---|
| 7 | import os, string, subprocess, sys |
---|
| 8 | import tempfile |
---|
| 9 | import re |
---|
| 10 | |
---|
| 11 | assert sys.version_info[:2] >= ( 2, 4 ) |
---|
| 12 | |
---|
| 13 | def stop_err(msg): |
---|
| 14 | sys.stderr.write(msg) |
---|
| 15 | sys.exit() |
---|
| 16 | |
---|
| 17 | def safe_bed_file(infile): |
---|
| 18 | """Make a BED file with track and browser lines ready for liftOver. |
---|
| 19 | |
---|
| 20 | liftOver will fail with track or browser lines. We can make it happy |
---|
| 21 | by converting these to comments. See: |
---|
| 22 | |
---|
| 23 | https://lists.soe.ucsc.edu/pipermail/genome/2007-May/013561.html |
---|
| 24 | """ |
---|
| 25 | fix_pat = re.compile("^(track|browser)") |
---|
| 26 | (fd, fname) = tempfile.mkstemp() |
---|
| 27 | in_handle = open(infile) |
---|
| 28 | out_handle = open(fname, "w") |
---|
| 29 | for line in in_handle: |
---|
| 30 | if fix_pat.match(line): |
---|
| 31 | line = "#" + line |
---|
| 32 | out_handle.write(line) |
---|
| 33 | in_handle.close() |
---|
| 34 | out_handle.close() |
---|
| 35 | return fname |
---|
| 36 | |
---|
| 37 | if len( sys.argv ) != 7: |
---|
| 38 | stop_err( "USAGE: prog input out_file1 out_file2 input_dbkey output_dbkey minMatch" ) |
---|
| 39 | |
---|
| 40 | infile = sys.argv[1] |
---|
| 41 | outfile1 = sys.argv[2] |
---|
| 42 | outfile2 = sys.argv[3] |
---|
| 43 | in_dbkey = sys.argv[4] |
---|
| 44 | mapfilepath = sys.argv[5] |
---|
| 45 | minMatch = sys.argv[6] |
---|
| 46 | try: |
---|
| 47 | assert float(minMatch) |
---|
| 48 | except: |
---|
| 49 | minMatch = 0.1 |
---|
| 50 | #ensure dbkey is set |
---|
| 51 | if in_dbkey == "?": |
---|
| 52 | stop_err( "Input dataset genome build unspecified, click the pencil icon in the history item to specify it." ) |
---|
| 53 | |
---|
| 54 | if not os.path.isfile( mapfilepath ): |
---|
| 55 | stop_err( "%s mapping is not currently available." % ( mapfilepath.split('/')[-1].split('.')[0] ) ) |
---|
| 56 | |
---|
| 57 | safe_infile = safe_bed_file(infile) |
---|
| 58 | cmd_line = "liftOver -minMatch=" + str(minMatch) + " " + safe_infile + " " + mapfilepath + " " + outfile1 + " " + outfile2 + " > /dev/null" |
---|
| 59 | try: |
---|
| 60 | # have to nest try-except in try-finally to handle 2.4 |
---|
| 61 | try: |
---|
| 62 | proc = subprocess.Popen( args=cmd_line, shell=True, stderr=subprocess.PIPE ) |
---|
| 63 | returncode = proc.wait() |
---|
| 64 | stderr = proc.stderr.read() |
---|
| 65 | if returncode != 0: |
---|
| 66 | raise Exception, stderr |
---|
| 67 | except Exception, e: |
---|
| 68 | raise Exception, 'Exception caught attempting conversion: ' + str( e ) |
---|
| 69 | finally: |
---|
| 70 | os.remove(safe_infile) |
---|