root/galaxy-central/tools/extract/liftOver_wrapper.py @ 2

リビジョン 2, 2.1 KB (コミッタ: hatakeyama, 14 年 前)

import galaxy-central

行番号 
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)
Note: リポジトリブラウザについてのヘルプは TracBrowser を参照してください。