annotate gff_fmap.py @ 5:6e589f267c14

Uploaded
author devteam
date Tue, 04 Nov 2014 12:15:19 -0500
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
5
6e589f267c14 Uploaded
devteam
parents:
diff changeset
1 #!/usr/bin/env python
6e589f267c14 Uploaded
devteam
parents:
diff changeset
2 """
6e589f267c14 Uploaded
devteam
parents:
diff changeset
3 GFF feature mapping program, to find the relation between different features described in a given GFF file.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
4
6e589f267c14 Uploaded
devteam
parents:
diff changeset
5 Usage:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
6 python gff_fmap.py in.gff > out.txt
6e589f267c14 Uploaded
devteam
parents:
diff changeset
7
6e589f267c14 Uploaded
devteam
parents:
diff changeset
8 Courtesy: Brad Chapman
6e589f267c14 Uploaded
devteam
parents:
diff changeset
9 Few functions are inherited from bcbio-GFFParser program.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
10 """
6e589f267c14 Uploaded
devteam
parents:
diff changeset
11
6e589f267c14 Uploaded
devteam
parents:
diff changeset
12 import re
6e589f267c14 Uploaded
devteam
parents:
diff changeset
13 import sys
6e589f267c14 Uploaded
devteam
parents:
diff changeset
14 import urllib
6e589f267c14 Uploaded
devteam
parents:
diff changeset
15 import collections
6e589f267c14 Uploaded
devteam
parents:
diff changeset
16 from helper import open_file
6e589f267c14 Uploaded
devteam
parents:
diff changeset
17
6e589f267c14 Uploaded
devteam
parents:
diff changeset
18 def _gff_line_map(line):
6e589f267c14 Uploaded
devteam
parents:
diff changeset
19 """Parses a line of GFF into a dictionary.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
20 Given an input line from a GFF file, this:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
21 - breaks it into component elements
6e589f267c14 Uploaded
devteam
parents:
diff changeset
22 - determines the type of attribute (flat, parent, child or annotation)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
23 - generates a dictionary of GFF info
6e589f267c14 Uploaded
devteam
parents:
diff changeset
24 """
6e589f267c14 Uploaded
devteam
parents:
diff changeset
25 gff3_kw_pat = re.compile("\w+=")
6e589f267c14 Uploaded
devteam
parents:
diff changeset
26 def _split_keyvals(keyval_str):
6e589f267c14 Uploaded
devteam
parents:
diff changeset
27 """Split key-value pairs in a GFF2, GTF and GFF3 compatible way.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
28 GFF3 has key value pairs like:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
29 count=9;gene=amx-2;sequence=SAGE:aacggagccg
6e589f267c14 Uploaded
devteam
parents:
diff changeset
30 GFF2 and GTF have:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
31 Sequence "Y74C9A" ; Note "Clone Y74C9A; Genbank AC024206"
6e589f267c14 Uploaded
devteam
parents:
diff changeset
32 name "fgenesh1_pg.C_chr_1000003"; transcriptId 869
6e589f267c14 Uploaded
devteam
parents:
diff changeset
33 """
6e589f267c14 Uploaded
devteam
parents:
diff changeset
34 quals = collections.defaultdict(list)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
35 if keyval_str is None:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
36 return quals
6e589f267c14 Uploaded
devteam
parents:
diff changeset
37 # ensembl GTF has a stray semi-colon at the end
6e589f267c14 Uploaded
devteam
parents:
diff changeset
38 if keyval_str[-1] == ';':
6e589f267c14 Uploaded
devteam
parents:
diff changeset
39 keyval_str = keyval_str[:-1]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
40 # GFF2/GTF has a semi-colon with at least one space after it.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
41 # It can have spaces on both sides; wormbase does this.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
42 # GFF3 works with no spaces.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
43 # Split at the first one we can recognize as working
6e589f267c14 Uploaded
devteam
parents:
diff changeset
44 parts = keyval_str.split(" ; ")
6e589f267c14 Uploaded
devteam
parents:
diff changeset
45 if len(parts) == 1:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
46 parts = keyval_str.split("; ")
6e589f267c14 Uploaded
devteam
parents:
diff changeset
47 if len(parts) == 1:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
48 parts = keyval_str.split(";")
6e589f267c14 Uploaded
devteam
parents:
diff changeset
49 # check if we have GFF3 style key-vals (with =)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
50 is_gff2 = True
6e589f267c14 Uploaded
devteam
parents:
diff changeset
51 if gff3_kw_pat.match(parts[0]):
6e589f267c14 Uploaded
devteam
parents:
diff changeset
52 is_gff2 = False
6e589f267c14 Uploaded
devteam
parents:
diff changeset
53 key_vals = [p.split('=') for p in parts]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
54 # otherwise, we are separated by a space with a key as the first item
6e589f267c14 Uploaded
devteam
parents:
diff changeset
55 else:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
56 pieces = []
6e589f267c14 Uploaded
devteam
parents:
diff changeset
57 for p in parts:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
58 # fix misplaced semi-colons in keys in some GFF2 files
6e589f267c14 Uploaded
devteam
parents:
diff changeset
59 if p and p[0] == ';':
6e589f267c14 Uploaded
devteam
parents:
diff changeset
60 p = p[1:]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
61 pieces.append(p.strip().split(" "))
6e589f267c14 Uploaded
devteam
parents:
diff changeset
62 key_vals = [(p[0], " ".join(p[1:])) for p in pieces]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
63 for key, val in key_vals:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
64 # remove quotes in GFF2 files
6e589f267c14 Uploaded
devteam
parents:
diff changeset
65 if (len(val) > 0 and val[0] == '"' and val[-1] == '"'):
6e589f267c14 Uploaded
devteam
parents:
diff changeset
66 val = val[1:-1]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
67 if val:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
68 quals[key].extend(val.split(','))
6e589f267c14 Uploaded
devteam
parents:
diff changeset
69 # if we don't have a value, make this a key=True/False style
6e589f267c14 Uploaded
devteam
parents:
diff changeset
70 # attribute
6e589f267c14 Uploaded
devteam
parents:
diff changeset
71 else:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
72 quals[key].append('true')
6e589f267c14 Uploaded
devteam
parents:
diff changeset
73 for key, vals in quals.items():
6e589f267c14 Uploaded
devteam
parents:
diff changeset
74 quals[key] = [urllib.unquote(v) for v in vals]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
75 return quals, is_gff2
6e589f267c14 Uploaded
devteam
parents:
diff changeset
76
6e589f267c14 Uploaded
devteam
parents:
diff changeset
77 def _nest_gff2_features(gff_parts):
6e589f267c14 Uploaded
devteam
parents:
diff changeset
78 """Provide nesting of GFF2 transcript parts with transcript IDs.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
79
6e589f267c14 Uploaded
devteam
parents:
diff changeset
80 exons and coding sequences are mapped to a parent with a transcript_id
6e589f267c14 Uploaded
devteam
parents:
diff changeset
81 in GFF2. This is implemented differently at different genome centers
6e589f267c14 Uploaded
devteam
parents:
diff changeset
82 and this function attempts to resolve that and map things to the GFF3
6e589f267c14 Uploaded
devteam
parents:
diff changeset
83 way of doing them.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
84 """
6e589f267c14 Uploaded
devteam
parents:
diff changeset
85 # map protein or transcript ids to a parent
6e589f267c14 Uploaded
devteam
parents:
diff changeset
86 for transcript_id in ["transcript_id", "transcriptId", "proteinId"]:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
87 try:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
88 gff_parts["quals"]["Parent"] = \
6e589f267c14 Uploaded
devteam
parents:
diff changeset
89 gff_parts["quals"][transcript_id]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
90 break
6e589f267c14 Uploaded
devteam
parents:
diff changeset
91 except KeyError:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
92 pass
6e589f267c14 Uploaded
devteam
parents:
diff changeset
93 # case for WormBase GFF -- everything labelled as Transcript or CDS
6e589f267c14 Uploaded
devteam
parents:
diff changeset
94 for flat_name in ["Transcript", "CDS"]:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
95 if gff_parts["quals"].has_key(flat_name):
6e589f267c14 Uploaded
devteam
parents:
diff changeset
96 # parent types
6e589f267c14 Uploaded
devteam
parents:
diff changeset
97 if gff_parts["type"] in [flat_name]:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
98 if not gff_parts["id"]:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
99 gff_parts["id"] = gff_parts["quals"][flat_name][0]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
100 gff_parts["quals"]["ID"] = [gff_parts["id"]]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
101 # children types
6e589f267c14 Uploaded
devteam
parents:
diff changeset
102 elif gff_parts["type"] in ["intron", "exon", "three_prime_UTR",
6e589f267c14 Uploaded
devteam
parents:
diff changeset
103 "coding_exon", "five_prime_UTR", "CDS", "stop_codon",
6e589f267c14 Uploaded
devteam
parents:
diff changeset
104 "start_codon"]:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
105 gff_parts["quals"]["Parent"] = gff_parts["quals"][flat_name]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
106 break
6e589f267c14 Uploaded
devteam
parents:
diff changeset
107
6e589f267c14 Uploaded
devteam
parents:
diff changeset
108 return gff_parts
6e589f267c14 Uploaded
devteam
parents:
diff changeset
109
6e589f267c14 Uploaded
devteam
parents:
diff changeset
110 line = line.strip()
6e589f267c14 Uploaded
devteam
parents:
diff changeset
111 if line == '':return [('directive', line)] # sometimes the blank lines will be there
6e589f267c14 Uploaded
devteam
parents:
diff changeset
112 if line[0] == '>':return [('directive', '')] # sometimes it will be a FATSA header
6e589f267c14 Uploaded
devteam
parents:
diff changeset
113 if line[0] == "#":
6e589f267c14 Uploaded
devteam
parents:
diff changeset
114 return [('directive', line[2:])]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
115 elif line:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
116 parts = line.split('\t')
6e589f267c14 Uploaded
devteam
parents:
diff changeset
117 if len(parts) == 1 and re.search(r'\w+', parts[0]):return [('directive', '')] ## GFF files with FASTA sequence together
6e589f267c14 Uploaded
devteam
parents:
diff changeset
118 assert len(parts) == 9, line
6e589f267c14 Uploaded
devteam
parents:
diff changeset
119 gff_parts = [(None if p == '.' else p) for p in parts]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
120 gff_info = dict()
6e589f267c14 Uploaded
devteam
parents:
diff changeset
121
6e589f267c14 Uploaded
devteam
parents:
diff changeset
122 # collect all of the base qualifiers for this item
6e589f267c14 Uploaded
devteam
parents:
diff changeset
123 quals, is_gff2 = _split_keyvals(gff_parts[8])
6e589f267c14 Uploaded
devteam
parents:
diff changeset
124
6e589f267c14 Uploaded
devteam
parents:
diff changeset
125 gff_info["is_gff2"] = is_gff2
6e589f267c14 Uploaded
devteam
parents:
diff changeset
126
6e589f267c14 Uploaded
devteam
parents:
diff changeset
127 if gff_parts[1]:quals["source"].append(gff_parts[1])
6e589f267c14 Uploaded
devteam
parents:
diff changeset
128 gff_info['quals'] = dict(quals)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
129
6e589f267c14 Uploaded
devteam
parents:
diff changeset
130 # if we are describing a location, then we are a feature
6e589f267c14 Uploaded
devteam
parents:
diff changeset
131 if gff_parts[3] and gff_parts[4]:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
132 gff_info['type'] = gff_parts[2]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
133 gff_info['id'] = quals.get('ID', [''])[0]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
134
6e589f267c14 Uploaded
devteam
parents:
diff changeset
135 if is_gff2:gff_info = _nest_gff2_features(gff_info)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
136 # features that have parents need to link so we can pick up
6e589f267c14 Uploaded
devteam
parents:
diff changeset
137 # the relationship
6e589f267c14 Uploaded
devteam
parents:
diff changeset
138 if gff_info['quals'].has_key('Parent'):
6e589f267c14 Uploaded
devteam
parents:
diff changeset
139 final_key = 'child'
6e589f267c14 Uploaded
devteam
parents:
diff changeset
140 elif gff_info['id']:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
141 final_key = 'parent'
6e589f267c14 Uploaded
devteam
parents:
diff changeset
142 # Handle flat features
6e589f267c14 Uploaded
devteam
parents:
diff changeset
143 else:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
144 final_key = 'feature'
6e589f267c14 Uploaded
devteam
parents:
diff changeset
145 # otherwise, associate these annotations with the full record
6e589f267c14 Uploaded
devteam
parents:
diff changeset
146 else:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
147 final_key = 'annotation'
6e589f267c14 Uploaded
devteam
parents:
diff changeset
148 return [(final_key, gff_info)]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
149
6e589f267c14 Uploaded
devteam
parents:
diff changeset
150 def parent_child_id_map(gff_handle):
6e589f267c14 Uploaded
devteam
parents:
diff changeset
151 """
6e589f267c14 Uploaded
devteam
parents:
diff changeset
152 Provide a mapping of parent to child relationships in the file.
6e589f267c14 Uploaded
devteam
parents:
diff changeset
153 Gives a dictionary of parent child relationships:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
154
6e589f267c14 Uploaded
devteam
parents:
diff changeset
155 keys -- tuple of (source, type) for each parent
6e589f267c14 Uploaded
devteam
parents:
diff changeset
156 values -- tuple of (source, type) as children of that parent
6e589f267c14 Uploaded
devteam
parents:
diff changeset
157 """
6e589f267c14 Uploaded
devteam
parents:
diff changeset
158 # collect all of the parent and child types mapped to IDs
6e589f267c14 Uploaded
devteam
parents:
diff changeset
159 parent_sts = dict()
6e589f267c14 Uploaded
devteam
parents:
diff changeset
160 child_sts = collections.defaultdict(list)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
161 for line in gff_handle:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
162 line_type, line_info = _gff_line_map(line)[0]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
163 if (line_type == 'parent' or (line_type == 'child' and line_info['id'])):
6e589f267c14 Uploaded
devteam
parents:
diff changeset
164 parent_sts[line_info['id']] = (line_info['quals']['source'][0], line_info['type'])
6e589f267c14 Uploaded
devteam
parents:
diff changeset
165 if line_type == 'child':
6e589f267c14 Uploaded
devteam
parents:
diff changeset
166 for parent_id in line_info['quals']['Parent']:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
167 child_sts[parent_id].append((line_info['quals']['source'][0], line_info['type']))
6e589f267c14 Uploaded
devteam
parents:
diff changeset
168 gff_handle.close()
6e589f267c14 Uploaded
devteam
parents:
diff changeset
169 # generate a dictionary of the unique final type relationships
6e589f267c14 Uploaded
devteam
parents:
diff changeset
170 pc_map = collections.defaultdict(list)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
171 for parent_id, parent_type in parent_sts.items():
6e589f267c14 Uploaded
devteam
parents:
diff changeset
172 for child_type in child_sts[parent_id]:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
173 pc_map[parent_type].append(child_type)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
174 pc_final_map = dict()
6e589f267c14 Uploaded
devteam
parents:
diff changeset
175 for ptype, ctypes in pc_map.items():
6e589f267c14 Uploaded
devteam
parents:
diff changeset
176 unique_ctypes = list(set(ctypes))
6e589f267c14 Uploaded
devteam
parents:
diff changeset
177 unique_ctypes.sort()
6e589f267c14 Uploaded
devteam
parents:
diff changeset
178 pc_final_map[ptype] = unique_ctypes
6e589f267c14 Uploaded
devteam
parents:
diff changeset
179 # some cases the GFF file represents a single feature type
6e589f267c14 Uploaded
devteam
parents:
diff changeset
180 if not pc_final_map:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
181 for fid, stypes in parent_sts.items():
6e589f267c14 Uploaded
devteam
parents:
diff changeset
182 pc_final_map[stypes] = dict()
6e589f267c14 Uploaded
devteam
parents:
diff changeset
183 # generate a report on feature id mapping in the file
6e589f267c14 Uploaded
devteam
parents:
diff changeset
184 print '+---------------------+---------------------------------+'
6e589f267c14 Uploaded
devteam
parents:
diff changeset
185 print '| Parent feature type | Associated child feature type(s)|'
6e589f267c14 Uploaded
devteam
parents:
diff changeset
186 print '+---------------------+---------------------------------+'
6e589f267c14 Uploaded
devteam
parents:
diff changeset
187 for key, value in pc_final_map.items():
6e589f267c14 Uploaded
devteam
parents:
diff changeset
188 print key[0], key[1]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
189 for child_to in value:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
190 print '\t\t\t|-',child_to[0], child_to[1]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
191 print '+---------------------+---------------------------------+'
6e589f267c14 Uploaded
devteam
parents:
diff changeset
192
6e589f267c14 Uploaded
devteam
parents:
diff changeset
193
6e589f267c14 Uploaded
devteam
parents:
diff changeset
194 if __name__=='__main__':
6e589f267c14 Uploaded
devteam
parents:
diff changeset
195
6e589f267c14 Uploaded
devteam
parents:
diff changeset
196 try:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
197 gff_file = sys.argv[1]
6e589f267c14 Uploaded
devteam
parents:
diff changeset
198 except:
6e589f267c14 Uploaded
devteam
parents:
diff changeset
199 print __doc__
6e589f267c14 Uploaded
devteam
parents:
diff changeset
200 sys.exit(-1)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
201
6e589f267c14 Uploaded
devteam
parents:
diff changeset
202 gff_handle = open_file(gff_file)
6e589f267c14 Uploaded
devteam
parents:
diff changeset
203 parent_child_id_map(gff_handle)