1 |
#!/usr/bin/env python |
2 |
# encoding: utf-8 |
3 |
""" |
4 |
contig_to_chr_coords.py |
5 |
|
6 |
Created by Cole Trapnell on 2008-09-05. |
7 |
Copyright (c) 2008 Cole Trapnell. All rights reserved. |
8 |
""" |
9 |
|
10 |
import sys |
11 |
import getopt |
12 |
|
13 |
|
14 |
help_message = ''' |
15 |
Takes the NCBI seq_contig file and maps contig coords to whole chromosome |
16 |
coords in a GTF, GFF, or BED file |
17 |
|
18 |
contig_to_chr_coords.py <format_flag> <seq_contig.md> <junctions.bed|transcripts.gff|transcripts.gtf> |
19 |
''' |
20 |
|
21 |
|
22 |
class Usage(Exception): |
23 |
def __init__(self, msg): |
24 |
self.msg = msg |
25 |
|
26 |
|
27 |
def main(argv=None): |
28 |
if argv is None: |
29 |
argv = sys.argv |
30 |
try: |
31 |
try: |
32 |
opts, args = getopt.getopt(argv[1:], "ho:vbg", ["help", "output=", "bed", "gff"]) |
33 |
except getopt.error, msg: |
34 |
raise Usage(msg) |
35 |
|
36 |
arg_is_splice = False |
37 |
arg_is_gff = False |
38 |
# option processing |
39 |
for option, value in opts: |
40 |
if option == "-v": |
41 |
verbose = True |
42 |
if option in ("-h", "--help"): |
43 |
raise Usage(help_message) |
44 |
if option in ("-o", "--output"): |
45 |
output = value |
46 |
if option in ("-b", "--bed"): |
47 |
arg_is_splice = True |
48 |
if option in ("-g", "--gff"): |
49 |
arg_is_gff = True |
50 |
|
51 |
if (arg_is_splice == False and arg_is_gff == False) or (arg_is_splice == True and arg_is_gff == True): |
52 |
print >> sys.stderr, "Error: please specify either -b or -g, but not both" |
53 |
raise Usage(help_message) |
54 |
|
55 |
if len(args) < 1: |
56 |
raise Usage(help_message) |
57 |
contig_to_chr_file = open(args[0]) |
58 |
|
59 |
contigs = {} |
60 |
|
61 |
for line in contig_to_chr_file.readlines(): |
62 |
if line[0] == "#": |
63 |
continue |
64 |
line = line.strip() |
65 |
cols = line.split('\t') |
66 |
if len(cols) < 9: |
67 |
continue |
68 |
chromosome = cols[1] |
69 |
group = cols[8] |
70 |
feature_name = cols[5] |
71 |
if not feature_name in ["start", "end"]: |
72 |
contigs[feature_name] = (chromosome, int(cols[2])) |
73 |
#print feature_name, chromosome, int(cols[2]) |
74 |
|
75 |
if arg_is_gff: |
76 |
gff_file = open(args[1]) |
77 |
|
78 |
lines = gff_file.readlines() |
79 |
print lines[0], |
80 |
for line in lines[1:]: |
81 |
line = line.strip() |
82 |
cols = line.split('\t') |
83 |
if len(cols) < 8: |
84 |
continue |
85 |
|
86 |
contig = cols[0] |
87 |
chr_fields = contig.split('|') |
88 |
refseq_id = chr_fields[3] |
89 |
contig = contigs.get(refseq_id) |
90 |
chr_name = contig[0] |
91 |
pipe_idx = chr_name.find('|') |
92 |
|
93 |
if pipe_idx != -1: |
94 |
chr_name = chr_name[:pipe_idx] |
95 |
if contig != None: |
96 |
#print line |
97 |
left_pos = contig[1] + int(cols[3]) |
98 |
right_pos = contig[1] + int(cols[4]) |
99 |
|
100 |
print "chr%s\tTopHat\tisland\t%d\t%d\t%s\t.\t.\t%s" % (chr_name, left_pos, right_pos, cols[5],cols[8]) |
101 |
#print >>sys.stderr, "%s\t%d\t%d\t%s\t%s\t%s\t%s" % (contig[0], left_pos, right_pos,cols[3],cols[6],cols[0],cols[1]) |
102 |
|
103 |
|
104 |
if arg_is_splice: |
105 |
splice_file = open(args[1]) |
106 |
|
107 |
lines = splice_file.readlines() |
108 |
print lines[0], |
109 |
for line in lines[1:]: |
110 |
line = line.strip() |
111 |
cols = line.split('\t') |
112 |
contig = cols[0] |
113 |
chr_fields = contig.split('|') |
114 |
refseq_id = chr_fields[3] |
115 |
contig = contigs.get(refseq_id) |
116 |
chr_name = contig[0] |
117 |
pipe_idx = chr_name.find('|') |
118 |
|
119 |
if pipe_idx != -1: |
120 |
chr_name = chr_name[:pipe_idx] |
121 |
if contig != None: |
122 |
#print line |
123 |
left_pos = contig[1] + int(cols[1]) |
124 |
right_pos = contig[1] + int(cols[2]) |
125 |
|
126 |
print "chr%s\t%d\t%d\t%s\t0\t%s\t%s\t%s\t255,0,0\t2\t1,1\t%s" % (chr_name, left_pos, right_pos, cols[3],cols[5],left_pos, right_pos,cols[11]) |
127 |
#print >>sys.stderr, "%s\t%d\t%d\t%s\t%s\t%s\t%s" % (contig[0], left_pos, right_pos,cols[3],cols[6],cols[0],cols[1]) |
128 |
except Usage, err: |
129 |
print >> sys.stderr, sys.argv[0].split("/")[-1] + ": " + str(err.msg) |
130 |
print >> sys.stderr, "\t for help use --help" |
131 |
return 2 |
132 |
|
133 |
|
134 |
if __name__ == "__main__": |
135 |
sys.exit(main()) |