-
Notifications
You must be signed in to change notification settings - Fork 0
/
gbk2gtf_simple.py
55 lines (47 loc) · 1.81 KB
/
gbk2gtf_simple.py
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
#!/usr/bin/env
from BCBio import GFF
from Bio.Seq import Seq
from Bio.SeqRecord import SeqRecord
from Bio.SeqFeature import SeqFeature, FeatureLocation
from Bio import SeqIO
import sys, re
# Functions
def grep(l,s):
return [i for i in l if s in i]
# Program
source = sys.argv[3]
frame = "."
feattype = "exon"
score = "."
feats=set()
exonnumber = 1
f = open(sys.argv[2],'w')
for rec in SeqIO.parse( sys.argv[1], "genbank"):
if rec.features:
for feature in rec.features:
if feature.type == "CDS":
feats.add(feature.type)
seqname = rec.id
if feature.strand == 1:
strand = "+"
elif feature.strand == -1:
strand = "-"
else:
strand = "."
start = feature.location.start.position+1
end = feature.location.end.position
if 'gene' in feature.qualifiers:
genename = feature.qualifiers['gene'][0]
transcriptname = genename + "-1"
if 'product' in feature.qualifiers:
annotation = feature.qualifiers['product'][0]
if 'db_xref' in feature.qualifiers:
db_xref= feature.qualifiers['db_xref']
geneid = grep(db_xref,'GeneID')[0].replace('GeneID:', '')
transcriptid = "T" + geneid
else:
continue
output = "%s\t%s\t%s\t%s\t%s\t%s\t%s\t%s\tgene_id \"%s\"; transcript_id \"%s\"; exon_number \"%s\"; gene_name \"%s\"; transcript_name \"%s\"\n" % (seqname, source, feattype, start, end, score, strand, frame, geneid, transcriptid, exonnumber, genename, transcriptname)
f.write(output)
f.close()
print ("Converted : " + rec.id)