gclib/tophat_cpp/contig_to_chr_coords.py

#!/usr/bin/env python
# encoding: utf-8
"""
contig_to_chr_coords.py

Created by Cole Trapnell on 2008-09-05.
Copyright (c) 2008 Cole Trapnell. All rights reserved.
"""

import sys
import getopt


help_message = '''
Takes the NCBI seq_contig file and maps contig coords to whole chromosome 
coords in a GTF, GFF, or BED file

    contig_to_chr_coords.py <format_flag> <seq_contig.md> <junctions.bed|transcripts.gff|transcripts.gtf>
'''


class Usage(Exception):
    def __init__(self, msg):
        self.msg = msg


def main(argv=None):
    if argv is None:
        argv = sys.argv
    try:
        try:
            opts, args = getopt.getopt(argv[1:], "ho:vbg", ["help", "output=", "bed", "gff"])
        except getopt.error, msg:
            raise Usage(msg)
    
        arg_is_splice = False
        arg_is_gff = False
        # option processing
        for option, value in opts:
            if option == "-v":
                verbose = True
            if option in ("-h", "--help"):
                raise Usage(help_message)
            if option in ("-o", "--output"):
                output = value
            if option in ("-b", "--bed"):
                arg_is_splice = True
            if option in ("-g", "--gff"):
                arg_is_gff = True
        
        if (arg_is_splice == False and arg_is_gff == False) or (arg_is_splice == True and arg_is_gff == True):
                        print >> sys.stderr, "Error: please specify either -b or -g, but not both"
                        raise Usage(help_message)
        
        if len(args) < 1:
            raise Usage(help_message)
        contig_to_chr_file = open(args[0])
        
        contigs = {}
        
        for line in contig_to_chr_file.readlines():
            if line[0] == "#":
                continue
            line = line.strip()
            cols = line.split('\t')
            if len(cols) < 9:
                continue
            chromosome = cols[1]
            group = cols[8]
            feature_name = cols[5]
            if not feature_name in ["start", "end"]:
                contigs[feature_name] = (chromosome, int(cols[2]))
                #print feature_name, chromosome, int(cols[2])
        
        if arg_is_gff:
            gff_file = open(args[1])
        
            lines = gff_file.readlines()
            print lines[0],
            for line in lines[1:]:
                line = line.strip()
                cols = line.split('\t')
                if len(cols) < 8:
                    continue
                    
                contig = cols[0]
                chr_fields = contig.split('|')
                refseq_id = chr_fields[3]
                contig = contigs.get(refseq_id)
                chr_name = contig[0]
                pipe_idx = chr_name.find('|')

                if pipe_idx != -1:
                    chr_name = chr_name[:pipe_idx]
                if contig != None:
                    #print line
                    left_pos = contig[1] + int(cols[3])
                    right_pos = contig[1] + int(cols[4])
                
                    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])
                    #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])
            
        
        if arg_is_splice:
            splice_file = open(args[1])
        
            lines = splice_file.readlines()
            print lines[0],
            for line in lines[1:]:
                line = line.strip()
                cols = line.split('\t')
                contig = cols[0]
                chr_fields = contig.split('|')
                refseq_id = chr_fields[3]
                contig = contigs.get(refseq_id)
                chr_name = contig[0]
                pipe_idx = chr_name.find('|')

                if pipe_idx != -1:
                    chr_name = chr_name[:pipe_idx]
                if contig != None:
                    #print line
                    left_pos = contig[1] + int(cols[1])
                    right_pos = contig[1] + int(cols[2])
                
                    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])
                    #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])
    except Usage, err:
        print >> sys.stderr, sys.argv[0].split("/")[-1] + ": " + str(err.msg)
        print >> sys.stderr, "\t for help use --help"
        return 2


if __name__ == "__main__":
    sys.exit(main())
Revision:	29
Committed:	Tue Aug 2 21:24:54 2011 UTC (13 years, 1 month ago) by gpertea
File size:	4631 byte(s)
Log Message:	adding tophat source work
Line	File contents
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())