· 8 years ago · Jan 24, 2018, 06:36 AM
1#!/bin/env python3
2
3# This file converts Cufflinks .gtf file to .bed file.
4import sys;
5import re;
6
7if len(sys.argv)<2:
8 print('This script converts Cufflinks predictions (.GTF) into .BED annotations.\n');
9 print('Usage: gtf2bed {OPTIONS} [.GTF file]\n');
10 print('Options:');
11 print('-c color\tSpecify the color of the track. This is a RGB value represented as "r,g,b". Default 255,0,0 (red)');
12 print('\nNote:');
13 print('1\tOnly "exon" and "transcript" are recognized in the feature field (3rd field).');
14 print('2\tIn the attribute list of .GTF file, the script tries to find "gene_id", "transcript_id" and "FPKM" attribute, and convert them as name and score field in .BED file.');
15 print('Author: Wei Li (li.david.wei AT gmail.com)');
16 sys.exit();
17
18color='255,0,0'
19
20
21for i in range(len(sys.argv)):
22 if sys.argv[i]=='-c':
23 color=sys.argv[i+1];
24
25
26
27def printbedline(estart,eend,field,nline):
28 try:
29 estp=estart[0]-1;
30 eedp=eend[-1];
31 # use regular expression to get transcript_id, gene_id and expression level
32 geneid=re.findall(r'gene_id \"([\w\.]+)\"',field[8])
33 transid=re.findall(r'transcript_id \"([\w\.]+)\"',field[8])
34 fpkmval=re.findall(r'FPKM \"([\d\.]+)\"',field[8])
35 if len(geneid)==0:
36 print('Warning: no gene_id field ',file=sys.stderr);
37 else:
38 geneid=geneid[0];
39 if len(transid)==0:
40 print('Warning: no transcript_id field',file=sys.stderr);
41 transid='Trans_'+str(nline);
42 else:
43 transid=transid[0];
44 if len(fpkmval)==0:
45 print('Warning: no FPKM field',file=sys.stderr);
46 fpkmval='100';
47 else:
48 fpkmval=fpkmval[0];
49 fpkmint=round(float(fpkmval));
50 print(field[0]+'\t'+str(estp)+'\t'+str(eedp)+'\t'+transid+'\t'+str(fpkmint)+'\t'+field[6]+'\t'+str(estp)+'\t'+str(eedp)+'\t'+color+'\t'+str(len(estart))+'\t',end='');
51 seglen=[eend[i]-estart[i]+1 for i in range(len(estart))];
52 segstart=[estart[i]-estart[0] for i in range(len(estart))];
53 strl=str(seglen[0]);
54 for i in range(1,len(seglen)):
55 strl+=','+str(seglen[i]);
56 strs=str(segstart[0]);
57 for i in range(1,len(segstart)):
58 strs+=','+str(segstart[i]);
59 print(strl+'\t'+strs);
60 except ValueError:
61 print('Error: non-number fields at line '+str(nline),file=sys.stderr);
62
63
64
65
66
67estart=[];
68eend=[];
69# read lines one to one
70nline=0;
71prevfield=[];
72for lines in open(sys.argv[-1]):
73 field=lines.strip().split('\t');
74 nline=nline+1;
75 if len(field)<9:
76 print('Error: the GTF should has at least 9 fields at line '+str(nline),file=sys.stderr);
77 continue;
78 if field[1]!='Cufflinks':
79 print('Warning: the second field is expected to be \'Cufflinks\' at line '+str(nline),file=sys.stderr);
80 if field[2]!='exon' and field[2] !='transcript':
81 print('Error: the third filed is expected to be \'exon\' or \'transcript\' at line '+str(nline),file=sys.stderr);
82 continue;
83 if field[2]=='exon':
84 try:
85 est=int(field[3]);
86 eed=int(field[4]);
87 estart+=[est];
88 eend+=[eed];
89 except ValueError:
90 print('Error: non-number fields at line '+str(nline),file=sys.stderr);
91 if field[2]=='transcript':
92 # A new transcript record, write
93 if len(estart)!=0:
94 printbedline(estart,eend,prevfield,nline);
95 prevfield=field;
96 estart=[];
97 eend=[];
98# the last record
99if len(estart)!=0:
100 printbedline(estart,eend,field,nline);