1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
|
"""
Readers extracting gene (exon and intron) information from bed / gtf / gff
formats.
- GeneReader: yields exons
- CDSReader: yields cds_exons
- FeatureReader: yields cds_exons, introns, exons
For gff/gtf, the start_codon stop_codon line types are merged with CDSs.
"""
import sys
from bx.bitset_utils import (
bitset_complement,
bitset_intersect,
bitset_subtract,
bitset_union,
)
def GeneReader(fh, format="gff"):
"""yield chrom, strand, gene_exons, name"""
known_formats = ("gff", "gtf", "bed")
if format not in known_formats:
print("{} format not in {}".format(format, ",".join(known_formats)), file=sys.stderr)
raise Exception("?")
if format == "bed":
for line in fh:
f = line.strip().split()
chrom = f[0]
chrom_start = int(f[1])
name = f[4]
strand = f[5]
int(f[6]) # cdsStart
int(f[7]) # cdsEnd
int(f[9]) # blockCount
blockSizes = [int(i) for i in f[10].strip(",").split(",")]
blockStarts = [chrom_start + int(i) for i in f[11].strip(",").split(",")]
# grab cdsStart - cdsEnd
gene_exons = []
for base, offset in zip(blockStarts, blockSizes):
exon_start = base
exon_end = base + offset
gene_exons.append((exon_start, exon_end))
yield chrom, strand, gene_exons, name
genelist = {}
grouplist = []
if format == "gff" or format == "gtf":
for line in fh:
if line.startswith("#"):
continue
fields = line.strip().split("\t")
if len(fields) < 9:
continue
# fields
chrom = fields[0]
ex_st = int(fields[3]) - 1 # make zero-centered
ex_end = int(fields[4]) # + 1 # make exclusive
strand = fields[6]
if format == "gtf":
group = fields[8].split(";")[0]
else:
group = fields[8]
if group not in grouplist:
grouplist.append(group)
if group not in genelist:
genelist[group] = (chrom, strand, [])
exons_i = 2
genelist[group][exons_i].append((ex_st, ex_end))
# for gene in genelist.values():
for gene in grouplist:
chrom, strand, gene_exons = genelist[gene]
gene_exons = bitset_union(gene_exons)
yield chrom, strand, gene_exons, gene
def CDSReader(fh, format="gff"):
"""yield chrom, strand, cds_exons, name"""
known_formats = ("gff", "gtf", "bed")
if format not in known_formats:
print("{} format not in {}".format(format, ",".join(known_formats)), file=sys.stderr)
raise Exception("?")
if format == "bed":
for line in fh:
f = line.strip().split()
chrom = f[0]
chrom_start = int(f[1])
name = f[4]
strand = f[5]
cdsStart = int(f[6])
cdsEnd = int(f[7])
int(f[9]) # blockCount
blockSizes = [int(i) for i in f[10].strip(",").split(",")]
blockStarts = [chrom_start + int(i) for i in f[11].strip(",").split(",")]
# grab cdsStart - cdsEnd
cds_exons = []
for base, offset in zip(blockStarts, blockSizes):
if (base + offset) < cdsStart:
continue
if base > cdsEnd:
continue
exon_start = max(base, cdsStart)
exon_end = min(base + offset, cdsEnd)
cds_exons.append((exon_start, exon_end))
yield chrom, strand, cds_exons, name
genelist = {}
grouplist = []
if format == "gff" or format == "gtf":
for line in fh:
if line.startswith("#"):
continue
fields = line.strip().split("\t")
if len(fields) < 9:
continue
if fields[2] not in ("CDS", "stop_codon", "start_codon"):
continue
# fields
chrom = fields[0]
ex_st = int(fields[3]) - 1 # make zero-centered
ex_end = int(fields[4]) # + 1 # make exclusive
strand = fields[6]
if format == "gtf":
group = fields[8].split(";")[0]
else:
group = fields[8]
if group not in grouplist:
grouplist.append(group)
if group not in genelist:
genelist[group] = (chrom, strand, [])
genelist[group][2].append((ex_st, ex_end))
# for gene in genelist.values():
for gene in grouplist:
chrom, strand, cds_exons = genelist[gene]
seqlen = sum(a[1] - a[0] for a in cds_exons)
overhang = seqlen % 3
if overhang > 0:
if strand == "+":
cds_exons[-1] = (cds_exons[-1][0], cds_exons[-1][1] - overhang)
else:
cds_exons[0] = (cds_exons[0][0] + overhang, cds_exons[0][1])
cds_exons = bitset_union(cds_exons)
yield chrom, strand, cds_exons, gene
def FeatureReader(fh, format="gff", alt_introns_subtract="exons", gtf_parse=None):
"""
yield chrom, strand, cds_exons, introns, exons, name
gtf_parse Example:
# parse gene_id from transcript_id "AC073130.2-001"; gene_id "TES";
gene_name = lambda s: s.split(';')[1].split()[1].strip('"')
for chrom, strand, cds_exons, introns, exons, name in FeatureReader( sys.stdin, format='gtf', gtf_parse=gene_name )
"""
known_formats = ("gff", "gtf", "bed")
if format not in known_formats:
print("{} format not in {}".format(format, ",".join(known_formats)), file=sys.stderr)
raise Exception("?")
if format == "bed":
for line in fh:
f = line.strip().split()
chrom = f[0]
chrom_start = int(f[1])
name = f[4]
strand = f[5]
cdsStart = int(f[6])
cdsEnd = int(f[7])
int(f[9]) # blockCount
blockSizes = [int(i) for i in f[10].strip(",").split(",")]
blockStarts = [chrom_start + int(i) for i in f[11].strip(",").split(",")]
# grab cdsStart - cdsEnd
cds_exons = []
exons = []
for base, offset in zip(blockStarts, blockSizes):
if (base + offset) < cdsStart:
continue
if base > cdsEnd:
continue
# exons
exon_start = base
exon_end = base + offset
exons.append((exon_start, exon_end))
# cds exons
exon_start = max(base, cdsStart)
exon_end = min(base + offset, cdsEnd)
cds_exons.append((exon_start, exon_end))
cds_exons = bitset_union(cds_exons)
exons = bitset_union(exons)
introns = bitset_complement(exons)
yield chrom, strand, cds_exons, introns, exons, name
genelist = {}
grouplist = []
if format == "gff" or format == "gtf":
for line in fh:
if line.startswith("#"):
continue
fields = line.strip().split("\t")
if len(fields) < 9:
continue
# fields
chrom = fields[0]
ex_st = int(fields[3]) - 1 # make zero-centered
ex_end = int(fields[4]) # + 1 # make exclusive
strand = fields[6]
if format == "gtf":
if not gtf_parse:
group = fields[8].split(";")[0]
else:
group = gtf_parse(fields[8])
else:
group = fields[8]
# Results are listed in the same order as encountered
if group not in grouplist:
grouplist.append(group)
if group not in genelist:
# chrom, strand, cds_exons, introns, exons, cds_start, cds_end
genelist[group] = [chrom, strand, [], [], [], None, None]
if fields[2] == "exon":
genelist[group][4].append((ex_st, ex_end))
elif fields[2] in ("CDS", "stop_codon", "start_codon"):
genelist[group][2].append((ex_st, ex_end))
if fields[2] == "start_codon":
if strand == "+":
genelist[group][5] = ex_st
else:
genelist[group][5] = ex_end
if fields[2] == "stop_codon":
if strand == "+":
genelist[group][5] = ex_end
else:
genelist[group][5] = ex_st
elif fields[2] == "intron":
genelist[group][3].append((ex_st, ex_end))
for gene in grouplist:
chrom, strand, cds_exons, introns, exons, cds_start, cds_end = genelist[gene]
cds_exons = bitset_union(cds_exons)
exons = bitset_union(exons)
# assure that cds exons were within the cds range
if cds_start is not None and cds_end is not None:
if strand == "+":
cds_exons = bitset_intersect(cds_exons, [(cds_start, cds_end)])
else:
cds_exons = bitset_intersect(cds_exons, [(cds_end, cds_start)])
# assure that introns are non-overlapping with themselves or exons
if alt_introns_subtract:
if alt_introns_subtract == "exons":
introns = bitset_subtract(introns, exons)
if alt_introns_subtract == "cds_exons":
introns = bitset_subtract(introns, cds_exons)
else:
introns = bitset_union(introns)
# assure CDS is a multiple of 3, trim from last exon if necessary
seqlen = sum(a[1] - a[0] for a in cds_exons)
overhang = seqlen % 3
if overhang > 0:
if strand == "+":
cds_exons[-1] = (cds_exons[-1][0], cds_exons[-1][1] - overhang)
else:
cds_exons[0] = (cds_exons[0][0] + overhang, cds_exons[0][1])
yield chrom, strand, cds_exons, introns, exons, gene
|
