PageRenderTime 24ms CodeModel.GetById 13ms app.highlight 7ms RepoModel.GetById 1ms app.codeStats 0ms

/tools/liftover/bedLiftover_wrapper.py

https://bitbucket.org/cistrome/cistrome-harvard/
Python | 70 lines | 48 code | 8 blank | 14 comment | 8 complexity | bc350105dddeea0035e33e9859e8a2bc MD5 | raw file
 1#!/usr/bin/env python
 2#Guruprasad Ananda
 3"""
 4Converts coordinates from one build/assembly to another using liftOver binary and mapping files downloaded from UCSC.
 5"""
 6
 7import os, string, subprocess, sys
 8import tempfile
 9import re
10
11assert sys.version_info[:2] >= ( 2, 4 )
12
13def stop_err(msg):
14    sys.stderr.write(msg)
15    sys.exit()
16
17def 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    
37if len( sys.argv ) != 7:
38    stop_err( "USAGE: prog input out_file1 out_file2 input_dbkey output_dbkey minMatch" )
39
40infile = sys.argv[1]
41outfile1 = sys.argv[2]
42outfile2 = sys.argv[3]
43in_dbkey = sys.argv[4]
44mapfilepath = sys.argv[5]
45minMatch = sys.argv[6]
46try:
47    assert float(minMatch)
48except:
49    minMatch = 0.1
50#ensure dbkey is set
51if in_dbkey == "?": 
52    stop_err( "Input dataset genome build unspecified, click the pencil icon in the history item to specify it." )
53
54if not os.path.isfile( mapfilepath ):
55    stop_err( "%s mapping is not currently available."  % ( mapfilepath.split('/')[-1].split('.')[0] ) )
56
57safe_infile = safe_bed_file(infile)
58cmd_line = "liftOver -minMatch=" + str(minMatch) + " " + safe_infile + " " + mapfilepath + " " + outfile1 + " " + outfile2 + "  > /dev/null"
59try:
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 )
69finally:
70    os.remove(safe_infile)