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) |
---|