1 |
gpertea |
29 |
#include "gff.h" |
2 |
|
|
|
3 |
|
|
//GffNames* GffReader::names=NULL; |
4 |
|
|
GffNames* GffObj::names=NULL; |
5 |
|
|
//global set of feature names, attribute names etc. |
6 |
|
|
// -- common for all GffObjs in current application! |
7 |
|
|
|
8 |
|
|
const uint GFF_MAX_LOCUS = 7000000; //longest known gene in human is ~2.2M, UCSC claims a gene for mouse of ~ 3.1 M |
9 |
|
|
const uint GFF_MAX_EXON = 30000; //longest known exon in human is ~11K |
10 |
|
|
const uint GFF_MAX_INTRON= 6000000; |
11 |
|
|
bool gff_show_warnings = false; //global setting, set by GffReader->showWarnings() |
12 |
|
|
const int gff_fid_mRNA=0; |
13 |
|
|
const int gff_fid_transcript=1; |
14 |
|
|
const int gff_fid_exon=2; |
15 |
|
|
const int gff_fid_CDS=3; //never really used in GffObj ftype_id or subftype_id |
16 |
|
|
const uint gfo_flag_HAS_ERRORS = 0x00000001; |
17 |
|
|
const uint gfo_flag_CHILDREN_PROMOTED= 0x00000002; |
18 |
|
|
const uint gfo_flag_IS_GENE = 0x00000004; |
19 |
|
|
const uint gfo_flag_IS_TRANSCRIPT = 0x00000008; |
20 |
|
|
const uint gfo_flag_FROM_GFF3 = 0x00000010; |
21 |
|
|
const uint gfo_flag_BY_EXON = 0x00000020; //created by subfeature (exon) directly |
22 |
|
|
const uint gfo_flag_DISCARDED = 0x00000100; |
23 |
|
|
const uint gfo_flag_LST_KEEP = 0x00000200; |
24 |
|
|
const uint gfo_flag_LEVEL_MSK = 0x00FF0000; |
25 |
|
|
const byte gfo_flagShift_LEVEL = 16; |
26 |
|
|
|
27 |
|
|
void gffnames_ref(GffNames* &n) { |
28 |
|
|
if (n==NULL) n=new GffNames(); |
29 |
|
|
n->numrefs++; |
30 |
|
|
} |
31 |
|
|
|
32 |
|
|
void gffnames_unref(GffNames* &n) { |
33 |
|
|
if (n==NULL) GError("Error: attempt to remove reference to null GffNames object!\n"); |
34 |
|
|
n->numrefs--; |
35 |
|
|
if (n->numrefs==0) { delete n; n=NULL; } |
36 |
|
|
} |
37 |
|
|
|
38 |
|
|
int gfo_cmpByLoc(const pointer p1, const pointer p2) { |
39 |
|
|
|
40 |
|
|
GffObj& g1=*((GffObj*)p1); |
41 |
|
|
GffObj& g2=*((GffObj*)p2); |
42 |
|
|
if (g1.gseq_id==g2.gseq_id) { |
43 |
|
|
if (g1.start!=g2.start) |
44 |
|
|
return (int)(g1.start-g2.start); |
45 |
|
|
else if (g1.getLevel()!=g2.getLevel()) |
46 |
|
|
return (int)(g1.getLevel()-g2.getLevel()); |
47 |
|
|
else |
48 |
|
|
if (g1.end!=g2.end) |
49 |
|
|
return (int)(g1.end-g2.end); |
50 |
|
|
else return strcmp(g1.getID(), g2.getID()); |
51 |
|
|
} |
52 |
|
|
else return (int)(g1.gseq_id-g2.gseq_id); |
53 |
|
|
} |
54 |
|
|
|
55 |
|
|
char* GffLine::extractAttr(const char* pre, bool caseStrict, bool enforce_GTF2) { |
56 |
|
|
//parse a key attribute and remove it from the info string |
57 |
|
|
//(only works for attributes that have values following them after ' ' or '=') |
58 |
|
|
static const char GTF2_ERR[]="Error parsing attribute %s ('\"' required) at GTF line:\n%s\n"; |
59 |
|
|
int lpre=strlen(pre); |
60 |
|
|
char cend=pre[lpre-1]; |
61 |
|
|
char* pos = (caseStrict) ? strstr(info, pre) : strifind(info, pre); |
62 |
|
|
if (pos==NULL) return NULL; |
63 |
|
|
char* findstart=info; |
64 |
|
|
//require word boundary on the left: |
65 |
|
|
while (pos!=NULL && pos!=info && *(pos-1)!=';' && *(pos-1)!=' ') { |
66 |
|
|
findstart=pos+lpre; |
67 |
|
|
pos = (caseStrict) ? strstr(findstart, pre) : strifind(findstart, pre); |
68 |
|
|
} |
69 |
|
|
if (pos==NULL) return NULL; |
70 |
|
|
if (cend!=' ' && cend!='=') { |
71 |
|
|
//require word boundary on the right: |
72 |
|
|
while (pos!=NULL && *(pos+lpre)!=' ' && *(pos+lpre)!='=') { |
73 |
|
|
findstart=pos+lpre; |
74 |
|
|
pos = (caseStrict) ? strstr(findstart, pre) : strifind(findstart, pre); |
75 |
|
|
} |
76 |
|
|
} |
77 |
|
|
if (pos==NULL) return NULL; |
78 |
|
|
char* vp=pos+lpre; |
79 |
|
|
while (*vp==' ') vp++; |
80 |
|
|
if (*vp==';' || *vp==0) |
81 |
|
|
GError("Error parsing value of GFF attribute \"%s\", line:\n%s\n", pre, dupline); |
82 |
|
|
bool dq_enclosed=false; //value string enclosed by double quotes |
83 |
|
|
if (*vp=='"') { |
84 |
|
|
dq_enclosed=true; |
85 |
|
|
vp++; |
86 |
|
|
} |
87 |
|
|
if (enforce_GTF2 && !dq_enclosed) |
88 |
|
|
GError(GTF2_ERR,pre, dupline); |
89 |
|
|
char* vend=vp; |
90 |
|
|
if (dq_enclosed) { |
91 |
|
|
while (*vend!='"' && *vend!=';' && *vend!=0) vend++; |
92 |
|
|
} |
93 |
|
|
else { |
94 |
|
|
while (*vend!=';' && *vend!=0) vend++; |
95 |
|
|
} |
96 |
|
|
if (enforce_GTF2 && *vend!='"') |
97 |
|
|
GError(GTF2_ERR, pre, dupline); |
98 |
|
|
char *r=Gstrdup(vp, vend-1); |
99 |
|
|
//-- now remove this attribute from the info string |
100 |
|
|
while (*vend!=0 && (*vend=='"' || *vend==';' || *vend==' ')) vend++; |
101 |
|
|
if (*vend==0) vend--; |
102 |
|
|
for (char *src=vend, *dest=pos;;src++,dest++) { |
103 |
|
|
*dest=*src; |
104 |
|
|
if (*src==0) break; |
105 |
|
|
} |
106 |
|
|
return r; |
107 |
|
|
} |
108 |
|
|
|
109 |
|
|
static char fnamelc[128]; |
110 |
|
|
|
111 |
|
|
GffLine::GffLine(GffReader* reader, const char* l) { |
112 |
|
|
llen=strlen(l); |
113 |
|
|
GMALLOC(line,llen+1); |
114 |
|
|
memcpy(line, l, llen+1); |
115 |
|
|
GMALLOC(dupline, llen+1); |
116 |
|
|
memcpy(dupline, l, llen+1); |
117 |
|
|
skip=true; |
118 |
|
|
gseqname=NULL; |
119 |
|
|
track=NULL; |
120 |
|
|
ftype=NULL; |
121 |
|
|
info=NULL; |
122 |
|
|
_parents=NULL; |
123 |
|
|
_parents_len=0; |
124 |
|
|
num_parents=0; |
125 |
|
|
parents=NULL; |
126 |
|
|
is_gff3=false; |
127 |
|
|
is_cds=false; |
128 |
|
|
is_transcript=false; |
129 |
|
|
is_exon=false; |
130 |
|
|
is_gene=false; |
131 |
|
|
exontype=0; |
132 |
|
|
gene_id=NULL; |
133 |
|
|
gene_name=NULL; |
134 |
|
|
qstart=0; |
135 |
|
|
qend=0; |
136 |
|
|
qlen=0; |
137 |
|
|
ID=NULL; |
138 |
|
|
char* t[9]; |
139 |
|
|
int i=0; |
140 |
|
|
int tidx=1; |
141 |
|
|
t[0]=line; |
142 |
|
|
|
143 |
|
|
while (line[i]!=0) { |
144 |
|
|
if (line[i]=='\t') { |
145 |
|
|
line[i]=0; |
146 |
|
|
t[tidx]=line+i+1; |
147 |
|
|
tidx++; |
148 |
|
|
if (tidx>8) break; |
149 |
|
|
} |
150 |
|
|
i++; |
151 |
|
|
} |
152 |
|
|
|
153 |
|
|
if (tidx<8) { // ignore non-GFF lines |
154 |
|
|
// GMessage("Warning: error parsing GFF/GTF line:\n%s\n", l); |
155 |
|
|
return; |
156 |
|
|
} |
157 |
|
|
gseqname=t[0]; |
158 |
|
|
track=t[1]; |
159 |
|
|
ftype=t[2]; |
160 |
|
|
info=t[8]; |
161 |
|
|
char* p=t[3]; |
162 |
|
|
if (!parseUInt(p,fstart)) |
163 |
|
|
GError("Error parsing start coordinate from GFF line:\n%s\n",l); |
164 |
|
|
p=t[4]; |
165 |
|
|
if (!parseUInt(p,fend)) |
166 |
|
|
GError("Error parsing end coordinate from GFF line:\n%s\n",l); |
167 |
|
|
if (fend<fstart) swap(fend,fstart); //make sure fstart>=fend, always |
168 |
|
|
p=t[5]; |
169 |
|
|
if (p[0]=='.' && p[1]==0) { |
170 |
|
|
score=0; |
171 |
|
|
} |
172 |
|
|
else { |
173 |
|
|
if (!parseDouble(p,score)) |
174 |
|
|
GError("Error parsing feature score from GFF line:\n%s\n",l); |
175 |
|
|
} |
176 |
|
|
strand=*t[6]; |
177 |
|
|
if (strand!='+' && strand!='-' && strand!='.') |
178 |
|
|
GError("Error parsing strand (%c) from GFF line:\n%s\n",strand,l); |
179 |
|
|
phase=*t[7]; // must be '.', '0', '1' or '2' |
180 |
|
|
ID=NULL; |
181 |
|
|
// exon/CDS/mrna filter |
182 |
|
|
strncpy(fnamelc, ftype, 127); |
183 |
|
|
fnamelc[127]=0; |
184 |
|
|
strlower(fnamelc); //convert to lower case |
185 |
|
|
bool is_t_data=false; |
186 |
|
|
if (strstr(fnamelc, "utr")!=NULL) { |
187 |
|
|
exontype=exgffUTR; |
188 |
|
|
is_exon=true; |
189 |
|
|
is_t_data=true; |
190 |
|
|
} |
191 |
|
|
else if (strstr(fnamelc, "exon")!=NULL) { |
192 |
|
|
exontype=exgffExon; |
193 |
|
|
is_exon=true; |
194 |
|
|
is_t_data=true; |
195 |
|
|
} |
196 |
|
|
else if (strstr(fnamelc, "stop") && |
197 |
|
|
(strstr(fnamelc, "codon") || strstr(fnamelc, "cds"))){ |
198 |
|
|
exontype=exgffStop; |
199 |
|
|
is_cds=true; //though some place it outside the last CDS segment |
200 |
|
|
is_t_data=true; |
201 |
|
|
} |
202 |
|
|
else if (strstr(fnamelc, "start") && |
203 |
|
|
((strstr(fnamelc, "codon")!=NULL) || strstr(fnamelc, "cds")!=NULL)){ |
204 |
|
|
exontype=exgffStart; |
205 |
|
|
is_cds=true; |
206 |
|
|
is_t_data=true; |
207 |
|
|
} |
208 |
|
|
else if (strcmp(fnamelc, "cds")==0) { |
209 |
|
|
exontype=exgffCDS; |
210 |
|
|
is_cds=true; |
211 |
|
|
is_t_data=true; |
212 |
|
|
} |
213 |
|
|
else if (endsWith(fnamelc, "gene") || startsWith(fnamelc, "gene")) { |
214 |
|
|
is_gene=true; |
215 |
|
|
is_t_data=true; //because its name will be attached to parented transcripts |
216 |
|
|
} |
217 |
|
|
else if (endsWith(fnamelc,"rna") || endsWith(fnamelc,"transcript")) { |
218 |
|
|
is_transcript=true; |
219 |
|
|
is_t_data=true; |
220 |
|
|
} |
221 |
|
|
|
222 |
|
|
if (reader->transcriptsOnly && !is_t_data) { |
223 |
|
|
char* id=extractAttr("ID="); |
224 |
|
|
if (id==NULL) id=extractAttr("transcript_id"); |
225 |
|
|
//GMessage("Discarding non-transcript line:\n%s\n",l); |
226 |
|
|
if (id!=NULL) { |
227 |
|
|
reader->discarded_ids.Add(id, new int(1)); |
228 |
|
|
GFREE(id); |
229 |
|
|
} |
230 |
|
|
return; //skip this line, unwanted feature name |
231 |
|
|
} |
232 |
|
|
ID=extractAttr("ID="); |
233 |
|
|
char* Parent=extractAttr("Parent="); |
234 |
|
|
is_gff3=(ID!=NULL || Parent!=NULL); |
235 |
|
|
if (is_gff3) { |
236 |
|
|
//parse as GFF3 |
237 |
|
|
if (ID!=NULL) { |
238 |
|
|
//has ID attr so it's likely to be a parent feature |
239 |
|
|
//look for explicit gene name |
240 |
|
|
gene_name=extractAttr("gene_name=",false); |
241 |
|
|
if (gene_name==NULL) { |
242 |
|
|
gene_name=extractAttr("geneName=",false); |
243 |
|
|
if (gene_name==NULL) { |
244 |
|
|
gene_name=extractAttr("gene_sym=",false); |
245 |
|
|
if (gene_name==NULL) { |
246 |
|
|
gene_name=extractAttr("gene=",false); |
247 |
|
|
} |
248 |
|
|
} |
249 |
|
|
} |
250 |
|
|
gene_id=extractAttr("geneID=",false); |
251 |
|
|
if (gene_id==NULL) { |
252 |
|
|
gene_id=extractAttr("gene_id=",false); |
253 |
|
|
} |
254 |
|
|
if (is_gene) { |
255 |
|
|
//special case: keep the Name and ID attributes of the gene feature |
256 |
|
|
if (gene_name==NULL) |
257 |
|
|
gene_name=extractAttr("Name="); |
258 |
|
|
if (gene_id==NULL) //the ID is also gene_id in this case |
259 |
|
|
gene_id=Gstrdup(ID); |
260 |
|
|
//skip=false; |
261 |
|
|
//return; |
262 |
|
|
GFREE(Parent); //TMI, we really don't care about gene Parents? |
263 |
|
|
} //gene feature |
264 |
|
|
}// has GFF3 ID |
265 |
|
|
if (Parent!=NULL) { |
266 |
|
|
//keep Parent attr |
267 |
|
|
//parse multiple parents |
268 |
|
|
num_parents=1; |
269 |
|
|
p=Parent; |
270 |
|
|
int last_delim_pos=-1; |
271 |
|
|
while (*p!=';' && *p!=0) { |
272 |
|
|
if (*p==',' && *(p+1)!=0 && *(p+1)!=';') { |
273 |
|
|
num_parents++; |
274 |
|
|
last_delim_pos=(p-Parent); |
275 |
|
|
} |
276 |
|
|
p++; |
277 |
|
|
} |
278 |
|
|
_parents_len=p-Parent+1; |
279 |
|
|
_parents=Parent; |
280 |
|
|
GMALLOC(parents, num_parents*sizeof(char*)); |
281 |
|
|
parents[0]=_parents; |
282 |
|
|
int i=1; |
283 |
|
|
if (last_delim_pos>0) { |
284 |
|
|
for (p=_parents+1;p<=_parents+last_delim_pos;p++) { |
285 |
|
|
if (*p==',') { |
286 |
|
|
char* ep=p-1; |
287 |
|
|
while (*ep==' ' && ep>_parents) ep--; |
288 |
|
|
*(ep+1)=0; //end the string there |
289 |
|
|
parents[i]=p+1; |
290 |
|
|
i++; |
291 |
|
|
} |
292 |
|
|
} |
293 |
|
|
} |
294 |
|
|
} //has Parent field |
295 |
|
|
} //GFF3 |
296 |
|
|
else { // GTF-like expected |
297 |
|
|
Parent=extractAttr("transcript_id"); |
298 |
|
|
if (Parent!=NULL) { //GTF2 format detected |
299 |
|
|
if (is_transcript) { |
300 |
|
|
// atypical GTF with a parent transcript line declared |
301 |
|
|
ID=Parent; |
302 |
|
|
Parent=NULL; |
303 |
|
|
} |
304 |
|
|
gene_id=extractAttr("gene_id"); // for GTF this is the only attribute accepted as geneID |
305 |
|
|
gene_name=extractAttr("gene_name"); |
306 |
|
|
if (gene_name==NULL) { |
307 |
|
|
gene_name=extractAttr("gene_sym"); |
308 |
|
|
if (gene_name==NULL) |
309 |
|
|
gene_name=extractAttr("gene"); |
310 |
|
|
} |
311 |
|
|
//prepare for parseAttr by adding '=' character instead of spaces for all attributes |
312 |
|
|
//after the attribute name |
313 |
|
|
p=info; |
314 |
|
|
bool noed=true; //not edited after the last delim |
315 |
|
|
bool nsp=false; //non-space found after last delim |
316 |
|
|
while (*p!=0) { |
317 |
|
|
if (*p==' ') { |
318 |
|
|
if (nsp && noed) { |
319 |
|
|
*p='='; |
320 |
|
|
noed=false; |
321 |
|
|
p++; |
322 |
|
|
continue; |
323 |
|
|
} |
324 |
|
|
} |
325 |
|
|
else nsp=true; //non-space |
326 |
|
|
if (*p==';') { noed=true; nsp=false; } |
327 |
|
|
p++; |
328 |
|
|
} |
329 |
|
|
} //GTF2 detected (no parent line) |
330 |
|
|
else {// Parent is NULL, check for jigsaw format or other pre-GTF2 format |
331 |
|
|
//char* fexon=strstr(fnamelc, "exon"); |
332 |
|
|
//if (fexon!=NULL) { |
333 |
|
|
if (exontype==exgffExon) { |
334 |
|
|
if (startsWith(track,"jigsaw")) { |
335 |
|
|
is_cds=true; |
336 |
|
|
strcpy(track,"jigsaw"); |
337 |
|
|
p=strchr(info,';'); |
338 |
|
|
if (p==NULL) { Parent=Gstrdup(info); info=NULL; } |
339 |
|
|
else { Parent=Gstrdup(info,p-1); |
340 |
|
|
info=p+1; |
341 |
|
|
} |
342 |
|
|
} |
343 |
|
|
} //exon feature? |
344 |
|
|
if (Parent==NULL && exontype>=exgffCDS && |
345 |
|
|
(i=strcspn(info,"; \t\n\r"))<=(int)(strlen(info)+1)) { |
346 |
|
|
//one word ID ? really desperate attempt to parse it here |
347 |
|
|
Parent=Gstrdup(info,info+i-1); |
348 |
|
|
info=NULL; //discard anything else on the line |
349 |
|
|
} |
350 |
|
|
} |
351 |
|
|
if (Parent!=NULL) { //GTF transcript_id for exon/CDS feature |
352 |
|
|
_parents=Parent; |
353 |
|
|
GMALLOC(parents,sizeof(char*)); |
354 |
|
|
num_parents=1; |
355 |
|
|
parents[0]=_parents; |
356 |
|
|
} |
357 |
|
|
} //GTF-like |
358 |
|
|
|
359 |
|
|
//parse other potentially useful features |
360 |
|
|
if (is_gff3) { |
361 |
|
|
if ((p=strstr(info,"Target="))!=NULL) { //has Target attr |
362 |
|
|
p+=7; |
363 |
|
|
while (*p!=';' && *p!=0 && *p!=' ') p++; |
364 |
|
|
if (*p!=' ') { |
365 |
|
|
GError("Error parsing target coordinates from GFF line:\n%s\n",l); |
366 |
|
|
} |
367 |
|
|
if (!parseUInt(p,qstart)) |
368 |
|
|
GError("Error parsing target start coordinate from GFF line:\n%s\n",l); |
369 |
|
|
if (*p!=' ') { |
370 |
|
|
GError("Error parsing next target coordinate from GFF line:\n%s\n",l); |
371 |
|
|
} |
372 |
|
|
p++; |
373 |
|
|
if (!parseUInt(p,qend)) |
374 |
|
|
GError("Error parsing target end coordinate from GFF line:\n%s\n",l); |
375 |
|
|
} |
376 |
|
|
if ((p=strifind(info,"Qreg="))!=NULL) { //has Qreg attr |
377 |
|
|
p+=5; |
378 |
|
|
if (!parseUInt(p,qstart)) |
379 |
|
|
GError("Error parsing target start coordinate from GFF line:\n%s\n",l); |
380 |
|
|
if (*p!='-') { |
381 |
|
|
GError("Error parsing next target coordinate from GFF line:\n%s\n",l); |
382 |
|
|
} |
383 |
|
|
p++; |
384 |
|
|
if (!parseUInt(p,qend)) |
385 |
|
|
GError("Error parsing target end coordinate from GFF line:\n%s\n",l); |
386 |
|
|
if (*p=='|' || *p==':') { |
387 |
|
|
p++; |
388 |
|
|
if (!parseUInt(p,qlen)) |
389 |
|
|
GError("Error parsing target length from GFF Qreg|: \n%s\n",l); |
390 |
|
|
} |
391 |
|
|
}//has Qreg attr |
392 |
|
|
if (qlen==0 && (p=strifind(info,"Qlen="))!=NULL) { |
393 |
|
|
p+=5; |
394 |
|
|
if (!parseUInt(p,qlen)) |
395 |
|
|
GError("Error parsing target length from GFF Qlen:\n%s\n",l); |
396 |
|
|
} |
397 |
|
|
}//parsing some useful attributes in GFF3 records |
398 |
|
|
if (ID==NULL && parents==NULL) { |
399 |
|
|
if (reader->gff_warns) |
400 |
|
|
GMessage("Warning: could not parse ID or Parent from GFF line:\n%s\n",dupline); |
401 |
|
|
return; //skip |
402 |
|
|
} |
403 |
|
|
skip=false; |
404 |
|
|
} |
405 |
|
|
|
406 |
|
|
int GffObj::addExon(GffReader* reader, GffLine* gl, bool keepAttr, bool noExonAttr) { |
407 |
|
|
//this will make sure we have the right subftype_id! |
408 |
|
|
int subf_id=-1; |
409 |
|
|
//if (ftype_id==gff_fid_mRNA) { //for mRNAs only parse known subfeatures! |
410 |
|
|
if (isTranscript()) { |
411 |
|
|
if (exon_ftype_id<0) {//exon_ftype_id=gff_fid_exon; |
412 |
|
|
if (gl->exontype>0) exon_ftype_id=gff_fid_exon; |
413 |
|
|
else exon_ftype_id=names->feats.addName(gl->ftype); |
414 |
|
|
} |
415 |
|
|
//any recognized mRNA segment gets the generic "exon" type (also applies to CDS) |
416 |
|
|
if (gl->exontype==0 && !gl->is_transcript) { |
417 |
|
|
//extraneous mRNA feature, discard |
418 |
|
|
if (reader->gff_warns) |
419 |
|
|
GMessage("Warning: discarding unrecognized transcript subfeature %s of %s\n", |
420 |
|
|
gl->ftype, gffID); |
421 |
|
|
return -1; |
422 |
|
|
} |
423 |
|
|
} |
424 |
|
|
else { //non-mRNA parent feature, check this subf type |
425 |
|
|
subf_id=names->feats.addName(gl->ftype); |
426 |
|
|
if (exon_ftype_id<0 || exons.Count()==0) //never assigned a subfeature type before (e.g. first exon being added) |
427 |
|
|
exon_ftype_id=subf_id; |
428 |
|
|
else { |
429 |
|
|
if (exon_ftype_id!=subf_id) { |
430 |
|
|
//if (subftype_id==ftype_id && exons.Count()==1 && exons[0]->start==start && exons[0]->end==end) { |
431 |
|
|
if (exon_ftype_id==ftype_id && exons.Count()==1 && exons[0]->start==start && exons[0]->end==end) { |
432 |
|
|
//the existing exon was just a dummy one created by default, discard it |
433 |
|
|
exons.Clear(); |
434 |
|
|
covlen=0; |
435 |
|
|
exon_ftype_id=subf_id; //allow the new subfeature to completely takeover |
436 |
|
|
} |
437 |
|
|
else { //multiple subfeatures, prefer those with |
438 |
|
|
if (reader->gff_warns) |
439 |
|
|
GMessage("GFF Warning: multiple subfeatures (%s and %s) found for %s, discarding ", |
440 |
|
|
names->feats.getName(subf_id), names->feats.getName(exon_ftype_id),gffID); |
441 |
|
|
if (gl->exontype!=0) { //new feature is an exon, discard previously parsed subfeatures |
442 |
|
|
if (reader->gff_warns) GMessage("%s.\n", names->feats.getName(exon_ftype_id)); |
443 |
|
|
exon_ftype_id=subf_id; |
444 |
|
|
exons.Clear(); |
445 |
|
|
covlen=0; |
446 |
|
|
} |
447 |
|
|
else { //discard new feature |
448 |
|
|
if (reader->gff_warns) GMessage("%s.\n", names->feats.getName(subf_id)); |
449 |
|
|
return -1; //skip this 2nd subfeature type for this parent! |
450 |
|
|
} |
451 |
|
|
} |
452 |
|
|
} //incoming subfeature is of different type |
453 |
|
|
} //new subfeature type |
454 |
|
|
} //non-mRNA parent |
455 |
|
|
int eidx=addExon(gl->fstart, gl->fend, gl->score, gl->phase, |
456 |
|
|
gl->qstart,gl->qend, gl->is_cds, gl->exontype); |
457 |
|
|
if (eidx<0) return eidx; //this should never happen |
458 |
|
|
if (keepAttr) { |
459 |
|
|
if (noExonAttr) { |
460 |
|
|
if (attrs==NULL) //place the parsed attributes directly at transcript level |
461 |
|
|
parseAttrs(attrs, gl->info); |
462 |
|
|
} |
463 |
|
|
else { //need all exon-level attributes |
464 |
|
|
parseAttrs(exons[eidx]->attrs, gl->info, true); |
465 |
|
|
} |
466 |
|
|
} |
467 |
|
|
return eidx; |
468 |
|
|
} |
469 |
|
|
|
470 |
|
|
|
471 |
|
|
int GffObj::addExon(uint segstart, uint segend, double sc, char fr, int qs, int qe, bool iscds, char exontype) { |
472 |
|
|
if (exons.Count()==0) { |
473 |
|
|
if (iscds) isCDS=true; //for now, assume CDS only if first "exon" given is a CDS |
474 |
|
|
if (exon_ftype_id<0) { |
475 |
|
|
exon_ftype_id = isTranscript() ? gff_fid_exon : ftype_id; |
476 |
|
|
} |
477 |
|
|
} |
478 |
|
|
//special treatment of start/stop codon features, they might be broken/split between exons |
479 |
|
|
//and in that case some providers will still give the wrong end coordinate as start+2 (e.g. UCSC) |
480 |
|
|
//so we should not trust the end coordinate for such features |
481 |
|
|
if (exontype==exgffStart || exontype==exgffStop) { |
482 |
|
|
if (strand=='-') segstart=segend; |
483 |
|
|
else segend=segstart; |
484 |
|
|
if (exontype==exgffStart) { |
485 |
|
|
if (CDstart==0 || segstart<CDstart) CDstart=segstart; |
486 |
|
|
} |
487 |
|
|
else { |
488 |
|
|
if (segstart>CDend) CDend=segstart; |
489 |
|
|
} |
490 |
|
|
} |
491 |
|
|
else if (iscds) { //update CDS anchors: |
492 |
|
|
if (CDstart==0 || segstart<CDstart) { |
493 |
|
|
CDstart=segstart; |
494 |
|
|
if (exontype==exgffCDS && strand=='+') CDphase=fr; |
495 |
|
|
} |
496 |
|
|
if (segend>CDend) { |
497 |
|
|
if (exontype==exgffCDS && strand=='-') CDphase=fr; |
498 |
|
|
CDend=segend; |
499 |
|
|
} |
500 |
|
|
} |
501 |
|
|
else { // not a CDS/start/stop |
502 |
|
|
isCDS=false; |
503 |
|
|
} |
504 |
|
|
if (qs || qe) { |
505 |
|
|
if (qs>qe) swap(qs,qe); |
506 |
|
|
if (qs==0) qs=1; |
507 |
|
|
} |
508 |
|
|
if (exontype>0) { //check for overlaps between exon-type segments |
509 |
|
|
int ovlen=0; |
510 |
|
|
int oi=exonOverlapIdx(segstart, segend, &ovlen); |
511 |
|
|
if (oi>=0) { //overlap existing segment |
512 |
|
|
if (ovlen==0) { |
513 |
|
|
//adjacent segments will be merged |
514 |
|
|
if ((exons[oi]->exontype==exgffUTR && exontype==exgffCDS) || |
515 |
|
|
(exons[oi]->exontype==exgffCDS && exontype==exgffUTR)) { |
516 |
|
|
expandExon(oi, segstart, segend, exgffCDSUTR, sc, fr, qs, qe); |
517 |
|
|
return oi; |
518 |
|
|
} |
519 |
|
|
} |
520 |
|
|
//only allow this for CDS within exon, stop_codon within exon, stop_codon within UTR, |
521 |
|
|
// start_codon within CDS or stop_codon within CDS |
522 |
|
|
if (exons[oi]->exontype>exontype && |
523 |
|
|
exons[oi]->start<=segstart && exons[oi]->end>=segend && |
524 |
|
|
!(exons[oi]->exontype==exgffUTR && exontype==exgffCDS)) { |
525 |
|
|
//larger segment given first, now the smaller included one is redundant |
526 |
|
|
return oi; //only used to store attributes from current GffLine |
527 |
|
|
} |
528 |
|
|
if (exontype>exons[oi]->exontype && |
529 |
|
|
segstart<=exons[oi]->start && segend>=exons[oi]->end && |
530 |
|
|
!(exontype==exgffUTR && exons[oi]->exontype==exgffCDS)) { |
531 |
|
|
//smaller segment given first, so we have to enlarge it |
532 |
|
|
expandExon(oi, segstart, segend, exontype, sc, fr, qs, qe); |
533 |
|
|
//this should also check for overlapping next exon (oi+1) ? |
534 |
|
|
return oi; |
535 |
|
|
} |
536 |
|
|
//there is also the special case of "ribosomal slippage exception" (programmed frameshift) |
537 |
|
|
//where two CDS segments may actually overlap for 1 or 2 bases, but there should be only one encompassing exon |
538 |
|
|
//if (ovlen>2 || exons[oi]->exontype!=exgffCDS || exontype!=exgffCDS) { |
539 |
|
|
// --> had to relax this because of some weird UCSC annotations with exons partially overlapping the CDS segments |
540 |
|
|
if (ovlen>2 && exons[oi]->exontype!=exgffUTR && exontype!=exgffUTR) { |
541 |
|
|
//important structural warning, will always print: |
542 |
|
|
if (gff_show_warnings) |
543 |
|
|
GMessage("GFF Warning: discarding overlapping feature segment (%d-%d) (vs %d-%d (%s)) for GFF ID %s on %s\n", |
544 |
|
|
segstart, segend, exons[oi]->start, exons[oi]->end, getSubfName(), gffID, getGSeqName()); |
545 |
|
|
hasErrors(true); |
546 |
|
|
return -1; //segment NOT added |
547 |
|
|
} |
548 |
|
|
// else add the segment if the overlap is small and between two CDS segments |
549 |
|
|
//TODO: we might want to add an attribute here with the slippage coordinate and size? |
550 |
|
|
}//overlap of existing segment |
551 |
|
|
} //check for overlap |
552 |
|
|
// --- no overlap, or accepted micro-overlap (ribosomal slippage) |
553 |
|
|
// create & add the new segment |
554 |
|
|
GffExon* enew=new GffExon(segstart, segend, sc, fr, qs, qe, exontype); |
555 |
|
|
int eidx=exons.Add(enew); |
556 |
|
|
if (eidx<0) { |
557 |
|
|
//this would actually be acceptable if the object is a "Gene" and "exons" are in fact isoforms |
558 |
|
|
if (gff_show_warnings) |
559 |
|
|
GMessage("GFF Warning: failed adding segment %d-%d for %s (discarded)!\n", |
560 |
|
|
segstart, segend, gffID); |
561 |
|
|
delete enew; |
562 |
|
|
hasErrors(true); |
563 |
|
|
return -1; |
564 |
|
|
} |
565 |
|
|
covlen+=(int)(exons[eidx]->end-exons[eidx]->start)+1; |
566 |
|
|
start=exons.First()->start; |
567 |
|
|
end=exons.Last()->end; |
568 |
|
|
if (uptr!=NULL) { //collect stats about the underlying genomic sequence |
569 |
|
|
GSeqStat* gsd=(GSeqStat*)uptr; |
570 |
|
|
if (start<gsd->mincoord) gsd->mincoord=start; |
571 |
|
|
if (end>gsd->maxcoord) gsd->maxcoord=end; |
572 |
|
|
if (this->len()>gsd->maxfeat_len) { |
573 |
|
|
gsd->maxfeat_len=this->len(); |
574 |
|
|
gsd->maxfeat=this; |
575 |
|
|
} |
576 |
|
|
} |
577 |
|
|
return eidx; |
578 |
|
|
} |
579 |
|
|
|
580 |
|
|
void GffObj::expandExon(int oi, uint segstart, uint segend, char exontype, double sc, char fr, int qs, int qe) { |
581 |
|
|
//oi is the index of the *first* overlapping segment found that must be enlarged |
582 |
|
|
covlen-=exons[oi]->len(); |
583 |
|
|
if (segstart<exons[oi]->start) |
584 |
|
|
exons[oi]->start=segstart; |
585 |
|
|
if (qs && qs<exons[oi]->qstart) exons[oi]->qstart=qs; |
586 |
|
|
if (segend>exons[oi]->end) |
587 |
|
|
exons[oi]->end=segend; |
588 |
|
|
if (qe && qe>exons[oi]->qend) exons[oi]->qend=qe; |
589 |
|
|
//warning: score cannot be properly adjusted! e.g. if it's a p-value it's just going to get worse |
590 |
|
|
if (sc!=0) exons[oi]->score=sc; |
591 |
|
|
covlen+=exons[oi]->len(); |
592 |
|
|
//if (exons[oi]->exontype< exontype) -- always true |
593 |
|
|
exons[oi]->exontype = exontype; |
594 |
|
|
if (exontype==exgffCDS) exons[oi]->phase=fr; |
595 |
|
|
//we must check if any more exons are also overlapping this |
596 |
|
|
int ni=oi+1; //next exon index after oi |
597 |
|
|
while (ni<exons.Count() && segend>=exons[ni]->start) { // next segment overlaps new enlarged segment |
598 |
|
|
//only allow this if next segment is fully included, and a subordinate |
599 |
|
|
if (exons[ni]->exontype<exontype && exons[ni]->end<=segend) { |
600 |
|
|
/* I guess we have to relax this due to stupid UCSC hg18 files having a start_codon sticking out |
601 |
|
|
chr1 hg18_knownGene start_codon 69806911 69806913 0.000000 + . |
602 |
|
|
chr1 hg18_knownGene CDS 69806911 69806912 0.000000 + 0 |
603 |
|
|
chr1 hg18_knownGene exon 69805456 69806912 0.000000 + . |
604 |
|
|
*/ |
605 |
|
|
if (exons[ni]->qstart<exons[oi]->qstart) exons[oi]->qstart=exons[ni]->qstart; |
606 |
|
|
if (exons[ni]->qend>exons[oi]->qend) exons[oi]->qend=exons[ni]->qend; |
607 |
|
|
exons.Delete(ni); |
608 |
|
|
} |
609 |
|
|
else { |
610 |
|
|
if (gff_show_warnings) GMessage("GFF Warning: overlapping existing exon(%d-%d) while expanding to %d-%d for GFF ID %s\n", |
611 |
|
|
exons[ni]->start, exons[ni]->end, segstart, segend, gffID); |
612 |
|
|
//hasErrors(true); |
613 |
|
|
break; |
614 |
|
|
} |
615 |
|
|
} |
616 |
|
|
// -- make sure any other related boundaries are updated: |
617 |
|
|
start=exons.First()->start; |
618 |
|
|
end=exons.Last()->end; |
619 |
|
|
if (uptr!=NULL) { //collect stats about the underlying genomic sequence |
620 |
|
|
GSeqStat* gsd=(GSeqStat*)uptr; |
621 |
|
|
if (start<gsd->mincoord) gsd->mincoord=start; |
622 |
|
|
if (end>gsd->maxcoord) gsd->maxcoord=end; |
623 |
|
|
if (this->len()>gsd->maxfeat_len) { |
624 |
|
|
gsd->maxfeat_len=this->len(); |
625 |
|
|
gsd->maxfeat=this; |
626 |
|
|
} |
627 |
|
|
} |
628 |
|
|
} |
629 |
|
|
|
630 |
|
|
void GffObj::removeExon(int idx) { |
631 |
|
|
/* |
632 |
|
|
if (idx==0 && segs[0].start==gstart) |
633 |
|
|
gstart=segs[1].start; |
634 |
|
|
if (idx==segcount && segs[segcount].end==gend) |
635 |
|
|
gend=segs[segcount-1].end; |
636 |
|
|
*/ |
637 |
|
|
if (idx<0 || idx>=exons.Count()) return; |
638 |
|
|
int segstart=exons[idx]->start; |
639 |
|
|
int segend=exons[idx]->end; |
640 |
|
|
exons.Delete(idx); |
641 |
|
|
covlen -= (int)(segend-segstart)+1; |
642 |
|
|
start=exons.First()->start; |
643 |
|
|
end=exons.Last()->end; |
644 |
|
|
if (isCDS) { CDstart=start; CDend=end; } |
645 |
|
|
} |
646 |
|
|
|
647 |
|
|
void GffObj::removeExon(GffExon* p) { |
648 |
|
|
for (int idx=0;idx<exons.Count();idx++) { |
649 |
|
|
if (exons[idx]==p) { |
650 |
|
|
int segstart=exons[idx]->start; |
651 |
|
|
int segend=exons[idx]->end; |
652 |
|
|
exons.Delete(idx); |
653 |
|
|
covlen -= (int)(segend-segstart)+1; |
654 |
|
|
start=exons.First()->start; |
655 |
|
|
end=exons.Last()->end; |
656 |
|
|
if (isCDS) { CDstart=start; CDend=end; } |
657 |
|
|
return; |
658 |
|
|
} |
659 |
|
|
} |
660 |
|
|
} |
661 |
|
|
|
662 |
|
|
|
663 |
|
|
|
664 |
|
|
GffObj::GffObj(GffReader *gfrd, GffLine* gffline, bool keepAttr, bool noExonAttr): |
665 |
|
|
GSeg(0,0), exons(true,true,false), children(1,false) { |
666 |
|
|
xstart=0; |
667 |
|
|
xend=0; |
668 |
|
|
xstatus=0; |
669 |
|
|
partial=false; |
670 |
|
|
isCDS=false; |
671 |
|
|
uptr=NULL; |
672 |
|
|
ulink=NULL; |
673 |
|
|
parent=NULL; |
674 |
|
|
udata=0; |
675 |
|
|
flags=0; |
676 |
|
|
CDstart=0; |
677 |
|
|
CDend=0; |
678 |
|
|
CDphase=0; |
679 |
|
|
geneID=NULL; |
680 |
|
|
gene_name=NULL; |
681 |
|
|
attrs=NULL; |
682 |
|
|
gffID=NULL; |
683 |
|
|
track_id=-1; |
684 |
|
|
gseq_id=-1; |
685 |
|
|
ftype_id=-1; |
686 |
|
|
exon_ftype_id=-1; |
687 |
|
|
strand='.'; |
688 |
|
|
if (gfrd==NULL) |
689 |
|
|
GError("Cannot use this GffObj constructor with a NULL GffReader!\n"); |
690 |
|
|
gffnames_ref(names); |
691 |
|
|
if (gfrd->names==NULL) gfrd->names=names; |
692 |
|
|
//qlen=0;qstart=0;qend=0; |
693 |
|
|
gscore=0; |
694 |
|
|
uscore=0; |
695 |
|
|
covlen=0; |
696 |
|
|
qcov=0; |
697 |
|
|
start=gffline->fstart; |
698 |
|
|
end=gffline->fend; |
699 |
|
|
gseq_id=names->gseqs.addName(gffline->gseqname); |
700 |
|
|
track_id=names->tracks.addName(gffline->track); |
701 |
|
|
strand=gffline->strand; |
702 |
|
|
qlen=gffline->qlen; |
703 |
|
|
qstart=gffline->qstart; |
704 |
|
|
qend=gffline->qend; |
705 |
|
|
//setup flags from gffline |
706 |
|
|
isCDS=gffline->is_cds; //for now |
707 |
|
|
isGene(gffline->is_gene); |
708 |
|
|
isTranscript(gffline->is_transcript || gffline->exontype!=0); |
709 |
|
|
fromGff3(gffline->is_gff3); |
710 |
|
|
|
711 |
|
|
if (gffline->parents!=NULL) { |
712 |
|
|
//GTF style -- create a GffObj directly by subfeature |
713 |
|
|
//(also possible orphan GFF3 exon line, or an exon given before its parent (chado)) |
714 |
|
|
if (gffline->exontype!=0) { //recognized exon-like feature |
715 |
|
|
ftype_id=gff_fid_transcript; //so this is some sort of transcript |
716 |
|
|
exon_ftype_id=gff_fid_exon; //subfeatures MUST be exons |
717 |
|
|
} |
718 |
|
|
else {//unrecognized subfeatures |
719 |
|
|
//make this GffObj of the same feature type |
720 |
|
|
ftype_id=names->feats.addName(gffline->ftype); |
721 |
|
|
} |
722 |
|
|
if (gffline->ID==NULL) { //typical GTF |
723 |
|
|
gffID=Gstrdup(gffline->parents[0]); |
724 |
|
|
this->createdByExon(true); |
725 |
|
|
//this is likely the first exon/segment of the feature |
726 |
|
|
addExon(gfrd, gffline, keepAttr, noExonAttr); |
727 |
|
|
} |
728 |
|
|
else { //a parented feature with an ID -- probably an orphan GFF3 line |
729 |
|
|
if (gffline->is_gff3 && gffline->exontype!=0) { |
730 |
|
|
//premature exon given before its parent transcript |
731 |
|
|
//create the transcript entry here |
732 |
|
|
gffID=Gstrdup(gffline->parents[0]); |
733 |
|
|
this->createdByExon(true); |
734 |
|
|
//this is the first exon/segment of the transcript |
735 |
|
|
addExon(gfrd, gffline, keepAttr, noExonAttr); |
736 |
|
|
} |
737 |
|
|
else { //unrecognized non-exon feature ? use the ID instead |
738 |
|
|
gffID=Gstrdup(gffline->ID); |
739 |
|
|
if (keepAttr) this->parseAttrs(attrs, gffline->info); |
740 |
|
|
} |
741 |
|
|
} |
742 |
|
|
} //subfeature given directly |
743 |
|
|
else { //gffline->parents==NULL |
744 |
|
|
//create a parent feature in its own right |
745 |
|
|
gscore=gffline->score; |
746 |
|
|
if (gffline->ID==NULL || gffline->ID[0]==0) |
747 |
|
|
GError("Error: no ID found for GFF record start\n"); |
748 |
|
|
gffID=Gstrdup(gffline->ID); //there must be an ID here |
749 |
|
|
//if (gffline->is_transcript) ftype_id=gff_fid_mRNA; |
750 |
|
|
//else |
751 |
|
|
ftype_id=names->feats.addName(gffline->ftype); |
752 |
|
|
if (gffline->is_transcript) |
753 |
|
|
exon_ftype_id=gff_fid_exon; |
754 |
|
|
|
755 |
|
|
if (keepAttr) this->parseAttrs(attrs, gffline->info); |
756 |
|
|
}//no parent |
757 |
|
|
|
758 |
|
|
if (gffline->gene_name!=NULL) { |
759 |
|
|
gene_name=Gstrdup(gffline->gene_name); |
760 |
|
|
} |
761 |
|
|
if (gffline->gene_id!=NULL) { |
762 |
|
|
geneID=Gstrdup(gffline->gene_id); |
763 |
|
|
} |
764 |
|
|
|
765 |
|
|
GSeqStat* gsd=gfrd->gseqstats.AddIfNew(new GSeqStat(gseq_id,names->gseqs.lastNameUsed()),true); |
766 |
|
|
uptr=gsd; |
767 |
|
|
if (start<gsd->mincoord) gsd->mincoord=start; |
768 |
|
|
if (end>gsd->maxcoord) gsd->maxcoord=end; |
769 |
|
|
if (this->len()>gsd->maxfeat_len) { |
770 |
|
|
gsd->maxfeat_len=this->len(); |
771 |
|
|
gsd->maxfeat=this; |
772 |
|
|
} |
773 |
|
|
} |
774 |
|
|
|
775 |
|
|
GffLine* GffReader::nextGffLine() { |
776 |
|
|
if (gffline!=NULL) return gffline; //caller should free gffline after processing |
777 |
|
|
while (gffline==NULL) { |
778 |
|
|
int llen=0; |
779 |
|
|
buflen=GFF_LINELEN-1; |
780 |
|
|
char* l=fgetline(linebuf, buflen, fh, &fpos, &llen); |
781 |
|
|
if (l==NULL) { |
782 |
|
|
return NULL; //end of file |
783 |
|
|
} |
784 |
|
|
int ns=0; //first nonspace position |
785 |
|
|
while (l[ns]!=0 && isspace(l[ns])) ns++; |
786 |
|
|
if (l[ns]=='#' || llen<10) continue; |
787 |
|
|
gffline=new GffLine(this, l); |
788 |
|
|
if (gffline->skip) { |
789 |
|
|
delete gffline; |
790 |
|
|
gffline=NULL; |
791 |
|
|
continue; |
792 |
|
|
} |
793 |
|
|
if (gffline->ID==NULL && gffline->parents==NULL) { //it must have an ID |
794 |
|
|
//this might not be needed, already checked in the GffLine constructor |
795 |
|
|
if (gff_warns) |
796 |
|
|
GMessage("Warning: malformed GFF line, no parent or record Id (kipping\n"); |
797 |
|
|
delete gffline; |
798 |
|
|
gffline=NULL; |
799 |
|
|
//continue; |
800 |
|
|
} |
801 |
|
|
} |
802 |
|
|
return gffline; |
803 |
|
|
} |
804 |
|
|
|
805 |
|
|
char* GffReader::gfoBuildId(const char* id, const char* ctg) { |
806 |
|
|
//caller must free the returned pointer |
807 |
|
|
char* buf=NULL; |
808 |
|
|
int idlen=strlen(id); |
809 |
|
|
GMALLOC(buf, idlen+strlen(ctg)+2); |
810 |
|
|
strcpy(buf, id); |
811 |
|
|
buf[idlen]='~'; |
812 |
|
|
strcpy(buf+idlen+1, ctg); |
813 |
|
|
return buf; |
814 |
|
|
} |
815 |
|
|
|
816 |
|
|
void GffReader::gfoRemove(const char* id, const char* ctg) { |
817 |
|
|
char* buf=gfoBuildId(id,ctg); |
818 |
|
|
phash.Remove(buf); |
819 |
|
|
GFREE(buf); |
820 |
|
|
} |
821 |
|
|
|
822 |
|
|
//Warning: if gflst gets altered, idx becomes obsolete |
823 |
|
|
GfoHolder* GffReader::gfoAdd(const char* id, const char* ctg, GffObj* gfo, int idx) { |
824 |
|
|
char* buf=gfoBuildId(id,ctg); |
825 |
|
|
GfoHolder* r=new GfoHolder(gfo,idx); |
826 |
|
|
phash.Add(buf, r); |
827 |
|
|
GFREE(buf); |
828 |
|
|
return r; |
829 |
|
|
} |
830 |
|
|
|
831 |
|
|
GfoHolder* GffReader::gfoFind(const char* id, const char* ctg) { |
832 |
|
|
char* buf=gfoBuildId(id,ctg); |
833 |
|
|
GfoHolder* r=phash.Find(buf); |
834 |
|
|
GFREE(buf); |
835 |
|
|
return r; |
836 |
|
|
} |
837 |
|
|
|
838 |
|
|
GfoHolder* GffReader::replaceGffRec(GffLine* gffline, bool keepAttr, bool noExonAttr, int replaceidx) { |
839 |
|
|
GffObj* newgfo=new GffObj(this, gffline, keepAttr, noExonAttr); |
840 |
|
|
GfoHolder* r=NULL; |
841 |
|
|
if (replaceidx>=0) { |
842 |
|
|
gflst.Put(replaceidx,newgfo); |
843 |
|
|
r=gfoAdd(newgfo->gffID, gffline->gseqname, newgfo, replaceidx); |
844 |
|
|
} |
845 |
|
|
else { |
846 |
|
|
int gfoidx=gflst.Add(newgfo); |
847 |
|
|
r=gfoAdd(newgfo->gffID, gffline->gseqname, newgfo, gfoidx); |
848 |
|
|
} |
849 |
|
|
if (gff_warns) { |
850 |
|
|
int* pcount=tids.Find(newgfo->gffID); |
851 |
|
|
if (pcount!=NULL) { |
852 |
|
|
if (gff_warns) GMessage("Warning: duplicate GFF ID: %s\n", newgfo->gffID); |
853 |
|
|
(*pcount)++; |
854 |
|
|
} |
855 |
|
|
else { |
856 |
|
|
tids.Add(newgfo->gffID,new int(1)); |
857 |
|
|
} |
858 |
|
|
} |
859 |
|
|
return r; |
860 |
|
|
} |
861 |
|
|
|
862 |
|
|
GfoHolder* GffReader::updateParent(GfoHolder* newgfh, GffObj* parent) { |
863 |
|
|
//assert(parent); |
864 |
|
|
//assert(newgfo); |
865 |
|
|
parent->children.Add(newgfh->gffobj); |
866 |
|
|
if (newgfh->gffobj->parent==NULL) newgfh->gffobj->parent=parent; |
867 |
|
|
newgfh->gffobj->setLevel(parent->getLevel()+1); |
868 |
|
|
if (parent->isGene()) { |
869 |
|
|
if (parent->gene_name!=NULL && newgfh->gffobj->gene_name==NULL) |
870 |
|
|
newgfh->gffobj->gene_name=Gstrdup(parent->gene_name); |
871 |
|
|
if (parent->geneID!=NULL && newgfh->gffobj->geneID==NULL) |
872 |
|
|
newgfh->gffobj->geneID=Gstrdup(parent->geneID); |
873 |
|
|
} |
874 |
|
|
|
875 |
|
|
return newgfh; |
876 |
|
|
} |
877 |
|
|
|
878 |
|
|
GfoHolder* GffReader::newGffRec(GffLine* gffline, bool keepAttr, bool noExonAttr, |
879 |
|
|
GffObj* parent, GffExon* pexon) { |
880 |
|
|
GffObj* newgfo=new GffObj(this, gffline, keepAttr, noExonAttr); |
881 |
|
|
GfoHolder* r=NULL; |
882 |
|
|
int gfoidx=gflst.Add(newgfo); |
883 |
|
|
r=gfoAdd(newgfo->gffID, gffline->gseqname, newgfo, gfoidx); |
884 |
|
|
if (parent!=NULL) { |
885 |
|
|
updateParent(r, parent); |
886 |
|
|
if (pexon!=NULL) parent->removeExon(pexon); |
887 |
|
|
} |
888 |
|
|
if (gff_warns) { |
889 |
|
|
int* pcount=tids.Find(newgfo->gffID); |
890 |
|
|
if (pcount!=NULL) { |
891 |
|
|
if (gff_warns) GMessage("Warning: duplicate GFF ID: %s\n", newgfo->gffID); |
892 |
|
|
(*pcount)++; |
893 |
|
|
} |
894 |
|
|
else { |
895 |
|
|
tids.Add(newgfo->gffID,new int(1)); |
896 |
|
|
} |
897 |
|
|
} |
898 |
|
|
return r; |
899 |
|
|
} |
900 |
|
|
|
901 |
|
|
GfoHolder* GffReader::updateGffRec(GfoHolder* prevgfo, GffLine* gffline, |
902 |
|
|
bool keepAttr) { |
903 |
|
|
if (prevgfo==NULL) return NULL; |
904 |
|
|
prevgfo->gffobj->createdByExon(false); |
905 |
|
|
prevgfo->gffobj->ftype_id=prevgfo->gffobj->names->feats.addName(gffline->ftype); |
906 |
|
|
prevgfo->gffobj->start=gffline->fstart; |
907 |
|
|
prevgfo->gffobj->end=gffline->fend; |
908 |
|
|
prevgfo->gffobj->isGene(gffline->is_gene); |
909 |
|
|
prevgfo->gffobj->isTranscript(gffline->is_transcript || gffline->exontype!=0); |
910 |
|
|
prevgfo->gffobj->fromGff3(gffline->is_gff3); |
911 |
|
|
if (keepAttr) { |
912 |
|
|
if (prevgfo->gffobj->attrs!=NULL) prevgfo->gffobj->attrs->Clear(); |
913 |
|
|
prevgfo->gffobj->parseAttrs(prevgfo->gffobj->attrs, gffline->info); |
914 |
|
|
} |
915 |
|
|
return prevgfo; |
916 |
|
|
} |
917 |
|
|
|
918 |
|
|
|
919 |
|
|
bool GffReader::addExonFeature(GfoHolder* prevgfo, GffLine* gffline, GHash<CNonExon>& pex, bool noExonAttr) { |
920 |
|
|
bool r=true; |
921 |
|
|
if (gffline->strand!=prevgfo->gffobj->strand) { |
922 |
|
|
GMessage("GFF Error: duplicate GFF ID '%s' (exons found on different strands of %s)\n", |
923 |
|
|
prevgfo->gffobj->gffID, prevgfo->gffobj->getGSeqName()); |
924 |
|
|
r=false; |
925 |
|
|
} |
926 |
|
|
int gdist=(gffline->fstart>prevgfo->gffobj->end) ? gffline->fstart-prevgfo->gffobj->end : |
927 |
|
|
((gffline->fend<prevgfo->gffobj->start)? prevgfo->gffobj->start-gffline->fend : |
928 |
|
|
0 ); |
929 |
|
|
if (gdist>(int)GFF_MAX_LOCUS) { //too far apart, most likely this is a duplicate ID |
930 |
|
|
GMessage("Error: duplicate GFF ID '%s' (or exons too far apart)!\n",prevgfo->gffobj->gffID); |
931 |
|
|
//validation_errors = true; |
932 |
|
|
r=false; |
933 |
|
|
if (!gff_warns) exit(1); |
934 |
|
|
} |
935 |
|
|
int eidx=prevgfo->gffobj->addExon(this, gffline, !noExonAttr, noExonAttr); |
936 |
|
|
if (eidx>=0 && gffline->ID!=NULL && gffline->exontype==0) |
937 |
|
|
subfPoolAdd(pex, prevgfo); |
938 |
|
|
return r; |
939 |
|
|
} |
940 |
|
|
|
941 |
|
|
CNonExon* GffReader::subfPoolCheck(GffLine* gffline, GHash<CNonExon>& pex, char*& subp_name) { |
942 |
|
|
CNonExon* subp=NULL; |
943 |
|
|
subp_name=NULL; |
944 |
|
|
for (int i=0;i<gffline->num_parents;i++) { |
945 |
|
|
if (transcriptsOnly && discarded_ids.Find(gffline->parents[i])!=NULL) |
946 |
|
|
continue; |
947 |
|
|
subp_name=gfoBuildId(gffline->parents[i], gffline->gseqname); //e.g. mRNA name |
948 |
|
|
subp=pex.Find(subp_name); |
949 |
|
|
if (subp!=NULL) |
950 |
|
|
return subp; |
951 |
|
|
GFREE(subp_name); |
952 |
|
|
} |
953 |
|
|
return NULL; |
954 |
|
|
} |
955 |
|
|
|
956 |
|
|
void GffReader::subfPoolAdd(GHash<CNonExon>& pex, GfoHolder* newgfo) { |
957 |
|
|
//this might become a parent feature later |
958 |
|
|
if (newgfo->gffobj->exons.Count()>0) { |
959 |
|
|
char* xbuf=gfoBuildId(gffline->ID, gffline->gseqname); |
960 |
|
|
pex.Add(xbuf, new CNonExon(newgfo->idx, newgfo->gffobj, |
961 |
|
|
newgfo->gffobj->exons[0], gffline)); |
962 |
|
|
GFREE(xbuf); |
963 |
|
|
} |
964 |
|
|
} |
965 |
|
|
|
966 |
|
|
GfoHolder* GffReader::promoteFeature(CNonExon* subp, char*& subp_name, GHash<CNonExon>& pex, |
967 |
|
|
bool keepAttr, bool noExonAttr) { |
968 |
|
|
GffObj* prevp=subp->parent; //grandparent of gffline (e.g. gene) |
969 |
|
|
if (prevp!=gflst[subp->idx]) |
970 |
|
|
GError("Error promoting subfeature %s, gflst index mismatch?!\n", subp->gffline->ID); |
971 |
|
|
subp->gffline->discardParent(); |
972 |
|
|
GfoHolder* gfoh=newGffRec(subp->gffline, keepAttr, noExonAttr, prevp, subp->exon); |
973 |
|
|
pex.Remove(subp_name); //no longer a potential parent, moved it to phash already |
974 |
|
|
prevp->promotedChildren(true); |
975 |
|
|
return gfoh; //returns the holder of newly promoted feature |
976 |
|
|
} |
977 |
|
|
|
978 |
|
|
//have to parse the whole file because exons can be scattered all over |
979 |
|
|
void GffReader::readAll(bool keepAttr, bool mergeCloseExons, bool noExonAttr) { |
980 |
|
|
bool validation_errors = false; |
981 |
|
|
//loc_debug=false; |
982 |
|
|
GHash<CNonExon> pex; //keep track of any "exon"-like features that have an ID |
983 |
|
|
//and thus could become promoted to parent features |
984 |
|
|
while (nextGffLine()!=NULL) { |
985 |
|
|
//seen this gff ID before? |
986 |
|
|
GfoHolder* prevseen=NULL; |
987 |
|
|
if (gffline->ID) //GFF3 |
988 |
|
|
prevseen=gfoFind(gffline->ID, gffline->gseqname); |
989 |
|
|
if (prevseen!=NULL) { |
990 |
|
|
if (prevseen->gffobj->createdByExon()) { |
991 |
|
|
updateGffRec(prevseen, gffline, keepAttr); |
992 |
|
|
} |
993 |
|
|
else { |
994 |
|
|
GMessage("Error: duplicate GFF ID '%s' encountered!\n",gffline->ID); |
995 |
|
|
validation_errors = true; |
996 |
|
|
if (gff_warns) { |
997 |
|
|
delete gffline; gffline=NULL; continue; |
998 |
|
|
} |
999 |
|
|
else exit(1); |
1000 |
|
|
} |
1001 |
|
|
} |
1002 |
|
|
if (gffline->parents==NULL) {//start GFF3-like record with no parent (mRNA, gene) |
1003 |
|
|
if (!prevseen) newGffRec(gffline, keepAttr, noExonAttr); |
1004 |
|
|
} |
1005 |
|
|
else { //--- it's a parented feature (could still be a mRNA) |
1006 |
|
|
bool found_parent=false; |
1007 |
|
|
GfoHolder* newgfo=prevseen; |
1008 |
|
|
for (int i=0;i<gffline->num_parents;i++) { |
1009 |
|
|
if (transcriptsOnly && discarded_ids.Find(gffline->parents[i])!=NULL) |
1010 |
|
|
continue; //skipping discarded parent feature |
1011 |
|
|
GfoHolder* parentgfo=gfoFind(gffline->parents[i], gffline->gseqname); |
1012 |
|
|
if (parentgfo!=NULL) { //parent GffObj parsed earlier |
1013 |
|
|
found_parent=true; |
1014 |
|
|
if (parentgfo->gffobj->isGene() && gffline->is_transcript |
1015 |
|
|
&& gffline->exontype==0) { |
1016 |
|
|
//not an exon, but a transcript parented by a gene |
1017 |
|
|
if (newgfo) { |
1018 |
|
|
updateParent(newgfo, parentgfo->gffobj); |
1019 |
|
|
} |
1020 |
|
|
else { |
1021 |
|
|
newgfo=newGffRec(gffline, keepAttr, noExonAttr, parentgfo->gffobj); |
1022 |
|
|
} |
1023 |
|
|
} |
1024 |
|
|
else { //potential exon subfeature |
1025 |
|
|
if (!addExonFeature(parentgfo, gffline, pex, noExonAttr)) |
1026 |
|
|
validation_errors=true; |
1027 |
|
|
} |
1028 |
|
|
} |
1029 |
|
|
} //for each parsed parent Id |
1030 |
|
|
if (!found_parent) { //new GTF-like record starting here with a subfeature directly |
1031 |
|
|
//or it could be some chado GFF3 barf with exons declared BEFORE their parent :( |
1032 |
|
|
//check if this feature isn't parented by a previously stored "exon" subfeature |
1033 |
|
|
char* subp_name=NULL; |
1034 |
|
|
CNonExon* subp=subfPoolCheck(gffline, pex, subp_name); |
1035 |
|
|
if (subp!=NULL) { //found a subfeature that is the parent of this gffline |
1036 |
|
|
//promote that subfeature to a full GffObj |
1037 |
|
|
GfoHolder* gfoh=promoteFeature(subp, subp_name, pex, keepAttr, noExonAttr); |
1038 |
|
|
//add current gffline as an exon of the newly promoted subfeature |
1039 |
|
|
if (!addExonFeature(gfoh, gffline, pex, noExonAttr)) |
1040 |
|
|
validation_errors=true; |
1041 |
|
|
} |
1042 |
|
|
else { //no parent seen before, create one directly with this exon |
1043 |
|
|
//loc_debug=true; |
1044 |
|
|
GfoHolder* newgfo=prevseen ? prevseen : newGffRec(gffline, keepAttr, noExonAttr); |
1045 |
|
|
if (gffline->ID!=NULL && gffline->exontype==0) |
1046 |
|
|
subfPoolAdd(pex, newgfo); |
1047 |
|
|
//even those with errors will be added here! |
1048 |
|
|
} |
1049 |
|
|
GFREE(subp_name); |
1050 |
|
|
} //no previous parent found |
1051 |
|
|
} //parented feature |
1052 |
|
|
//-- |
1053 |
|
|
delete gffline; |
1054 |
|
|
gffline=NULL; |
1055 |
|
|
}//while gff lines |
1056 |
|
|
gflst.finalize(this, mergeCloseExons, keepAttr, noExonAttr); //force sorting by locus if so constructed |
1057 |
|
|
// all gff records are now loaded in GList gflst |
1058 |
|
|
// so we can free the hash |
1059 |
|
|
phash.Clear(); |
1060 |
|
|
tids.Clear(); |
1061 |
|
|
if (validation_errors) { |
1062 |
|
|
exit(1); |
1063 |
|
|
} |
1064 |
|
|
} |
1065 |
|
|
|
1066 |
|
|
GffObj* GffObj::finalize(GffReader* gfr, bool mergeCloseExons, bool keepAttrs, bool noExonAttr) { |
1067 |
|
|
//merge |
1068 |
|
|
//always merge adjacent or overlapping segments |
1069 |
|
|
//but if mergeCloseExons then merge even when distance is up to 5 bases |
1070 |
|
|
udata=0; |
1071 |
|
|
uptr=NULL; |
1072 |
|
|
if (gfr->transcriptsOnly && !(isTranscript() || (isGene() && children.Count()==0))) { |
1073 |
|
|
isDiscarded(true); |
1074 |
|
|
} |
1075 |
|
|
if (ftype_id==gff_fid_transcript && CDstart>0) { |
1076 |
|
|
ftype_id=gff_fid_mRNA; |
1077 |
|
|
//exon_ftype_id=gff_fid_exon; |
1078 |
|
|
} |
1079 |
|
|
//if (ftype_id==gff_fid_mRNA || exon_ftype_id==gff_fid_exon || mergeCloseExons) { |
1080 |
|
|
if (isTranscript() || exon_ftype_id==gff_fid_exon || mergeCloseExons) { |
1081 |
|
|
int mindist=mergeCloseExons ? 5:1; |
1082 |
|
|
for (int i=0;i<exons.Count()-1;i++) { |
1083 |
|
|
int ni=i+1; |
1084 |
|
|
uint mend=exons[i]->end; |
1085 |
|
|
while (ni<exons.Count()) { |
1086 |
|
|
int dist=(int)(exons[ni]->start-mend); |
1087 |
|
|
if (dist>mindist) break; //no merging with next segment |
1088 |
|
|
if (gfr!=NULL && gfr->gff_warns && dist!=0 && (exons[ni]->exontype!=exgffUTR && exons[i]->exontype!=exgffUTR)) { |
1089 |
|
|
GMessage("GFF warning: merging adjacent/overlapping segments of %s on %s (%d-%d, %d-%d)\n", |
1090 |
|
|
gffID, getGSeqName(), exons[i]->start, exons[i]->end,exons[ni]->start, exons[ni]->end); |
1091 |
|
|
} |
1092 |
|
|
mend=exons[ni]->end; |
1093 |
|
|
covlen-=exons[i]->len(); |
1094 |
|
|
exons[i]->end=mend; |
1095 |
|
|
covlen+=exons[i]->len(); |
1096 |
|
|
covlen-=exons[ni]->len(); |
1097 |
|
|
if (exons[ni]->attrs!=NULL && (exons[i]->attrs==NULL || |
1098 |
|
|
exons[i]->attrs->Count()<exons[ni]->attrs->Count())) { |
1099 |
|
|
//use the other exon attributes, if more |
1100 |
|
|
delete(exons[i]->attrs); |
1101 |
|
|
exons[i]->attrs=exons[ni]->attrs; |
1102 |
|
|
exons[ni]->attrs=NULL; |
1103 |
|
|
} |
1104 |
|
|
exons.Delete(ni); |
1105 |
|
|
} //check for merge with next exon |
1106 |
|
|
} //for each exon |
1107 |
|
|
} |
1108 |
|
|
//attribute reduction for GTF records |
1109 |
|
|
if (keepAttrs && !noExonAttr && !fromGff3() |
1110 |
|
|
&& exons.Count()>0 && exons[0]->attrs!=NULL) { |
1111 |
|
|
bool attrs_discarded=false; |
1112 |
|
|
for (int a=0;a<exons[0]->attrs->Count();a++) { |
1113 |
|
|
int attr_name_id=exons[0]->attrs->Get(a)->attr_id; |
1114 |
|
|
char* attr_name=names->attrs.getName(attr_name_id); |
1115 |
|
|
char* attr_val =exons[0]->attrs->Get(a)->attr_val; |
1116 |
|
|
bool sameExonAttr=true; |
1117 |
|
|
for (int i=1;i<exons.Count();i++) { |
1118 |
|
|
char* ov=exons[i]->getAttr(attr_name_id); |
1119 |
|
|
if (ov==NULL || (strcmp(ov,attr_val)!=0)) { |
1120 |
|
|
sameExonAttr=false; |
1121 |
|
|
break; |
1122 |
|
|
} |
1123 |
|
|
} |
1124 |
|
|
if (sameExonAttr) { |
1125 |
|
|
//delete this attribute from exons level |
1126 |
|
|
attrs_discarded=true; |
1127 |
|
|
this->addAttr(attr_name, attr_val); |
1128 |
|
|
for (int i=1;i<exons.Count();i++) { |
1129 |
|
|
removeExonAttr(*(exons[i]), attr_name_id); |
1130 |
|
|
} |
1131 |
|
|
exons[0]->attrs->freeItem(a); |
1132 |
|
|
} |
1133 |
|
|
} |
1134 |
|
|
if (attrs_discarded) exons[0]->attrs->Pack(); |
1135 |
|
|
} |
1136 |
|
|
return this; |
1137 |
|
|
} |
1138 |
|
|
|
1139 |
|
|
void GffObj::parseAttrs(GffAttrs*& atrlist, char* info, bool isExon) { |
1140 |
|
|
if (names==NULL) |
1141 |
|
|
GError(ERR_NULL_GFNAMES, "parseAttrs()"); |
1142 |
|
|
if (atrlist==NULL) |
1143 |
|
|
atrlist=new GffAttrs(); |
1144 |
|
|
char* endinfo=info+strlen(info); |
1145 |
|
|
char* start=info; |
1146 |
|
|
char* pch=start; |
1147 |
|
|
while (start<endinfo) { |
1148 |
|
|
//skip spaces |
1149 |
|
|
while (*start==' ' && start<endinfo) start++; |
1150 |
|
|
pch=strchr(start, ';'); |
1151 |
|
|
if (pch==NULL) pch=endinfo; |
1152 |
|
|
else { |
1153 |
|
|
*pch='\0'; |
1154 |
|
|
pch++; |
1155 |
|
|
} |
1156 |
|
|
char* ech=strchr(start,'='); |
1157 |
|
|
if (ech!=NULL) { // attr=value format found |
1158 |
|
|
*ech='\0'; |
1159 |
|
|
//if (noExonAttr && (strcmp(start, "exon_number")==0 || strcmp(start, "exon")==0)) { start=pch; continue; } |
1160 |
|
|
if (strcmp(start, "exon_number")==0 || strcmp(start, "exon")==0 || |
1161 |
|
|
strcmp(start, "exon_id")==0) |
1162 |
|
|
{ start=pch; continue; } |
1163 |
|
|
ech++; |
1164 |
|
|
while (*ech==' ' && ech<endinfo) ech++;//skip extra spaces after the '=' |
1165 |
|
|
//atrlist->Add(new GffAttr(names->attrs.addName(start),ech)); |
1166 |
|
|
//make sure we don't add the same attribute more than once |
1167 |
|
|
if (isExon && (strcmp(start, "protein_id")==0)) { |
1168 |
|
|
//Ensembl special case |
1169 |
|
|
this->addAttr(start, ech); |
1170 |
|
|
start=pch; |
1171 |
|
|
continue; |
1172 |
|
|
} |
1173 |
|
|
atrlist->add_or_update(names, start, ech); |
1174 |
|
|
} |
1175 |
|
|
/* |
1176 |
|
|
else { //not an attr=value format |
1177 |
|
|
atrlist->Add(new GffAttr(names->attrs.addName(start),"1")); |
1178 |
|
|
} |
1179 |
|
|
*/ |
1180 |
|
|
start=pch; |
1181 |
|
|
} |
1182 |
|
|
if (atrlist->Count()==0) { delete atrlist; atrlist=NULL; } |
1183 |
|
|
} |
1184 |
|
|
|
1185 |
|
|
void GffObj::addAttr(const char* attrname, const char* attrvalue) { |
1186 |
|
|
if (this->attrs==NULL) |
1187 |
|
|
this->attrs=new GffAttrs(); |
1188 |
|
|
//this->attrs->Add(new GffAttr(names->attrs.addName(attrname),attrvalue)); |
1189 |
|
|
this->attrs->add_or_update(names, attrname, attrvalue); |
1190 |
|
|
} |
1191 |
|
|
|
1192 |
|
|
|
1193 |
|
|
void GffObj::setFeatureName(const char* feature) { |
1194 |
|
|
//change the feature name/type for a transcript |
1195 |
|
|
int fid=names->feats.addName(feature); |
1196 |
|
|
if (monoFeature() && exons.Count()>0) |
1197 |
|
|
this->exon_ftype_id=fid; |
1198 |
|
|
this->ftype_id=fid; |
1199 |
|
|
} |
1200 |
|
|
|
1201 |
|
|
void GffObj::setRefName(const char* newname) { |
1202 |
|
|
//change the feature name/type for a transcript |
1203 |
|
|
int rid=names->gseqs.addName(newname); |
1204 |
|
|
this->gseq_id=rid; |
1205 |
|
|
} |
1206 |
|
|
|
1207 |
|
|
|
1208 |
|
|
|
1209 |
|
|
int GffObj::removeAttr(const char* attrname, const char* attrval) { |
1210 |
|
|
if (this->attrs==NULL || attrname==NULL || attrname[0]==0) return 0; |
1211 |
|
|
int aid=this->names->attrs.getId(attrname); |
1212 |
|
|
if (aid<0) return 0; |
1213 |
|
|
int delcount=0; //could be more than one ? |
1214 |
|
|
for (int i=0;i<this->attrs->Count();i++) { |
1215 |
|
|
if (aid==this->attrs->Get(i)->attr_id) { |
1216 |
|
|
if (attrval==NULL || |
1217 |
|
|
strcmp(attrval, this->attrs->Get(i)->attr_val)==0) { |
1218 |
|
|
delcount++; |
1219 |
|
|
this->attrs->freeItem(i); |
1220 |
|
|
} |
1221 |
|
|
} |
1222 |
|
|
} |
1223 |
|
|
if (delcount>0) this->attrs->Pack(); |
1224 |
|
|
return delcount; |
1225 |
|
|
} |
1226 |
|
|
|
1227 |
|
|
int GffObj::removeAttr(int aid, const char* attrval) { |
1228 |
|
|
if (this->attrs==NULL || aid<0) return 0; |
1229 |
|
|
int delcount=0; //could be more than one ? |
1230 |
|
|
for (int i=0;i<this->attrs->Count();i++) { |
1231 |
|
|
if (aid==this->attrs->Get(i)->attr_id) { |
1232 |
|
|
if (attrval==NULL || |
1233 |
|
|
strcmp(attrval, this->attrs->Get(i)->attr_val)==0) { |
1234 |
|
|
delcount++; |
1235 |
|
|
this->attrs->freeItem(i); |
1236 |
|
|
} |
1237 |
|
|
} |
1238 |
|
|
} |
1239 |
|
|
if (delcount>0) this->attrs->Pack(); |
1240 |
|
|
return delcount; |
1241 |
|
|
} |
1242 |
|
|
|
1243 |
|
|
|
1244 |
|
|
int GffObj::removeExonAttr(GffExon& exon, const char* attrname, const char* attrval) { |
1245 |
|
|
if (exon.attrs==NULL || attrname==NULL || attrname[0]==0) return 0; |
1246 |
|
|
int aid=this->names->attrs.getId(attrname); |
1247 |
|
|
if (aid<0) return 0; |
1248 |
|
|
int delcount=0; //could be more than one |
1249 |
|
|
for (int i=0;i<exon.attrs->Count();i++) { |
1250 |
|
|
if (aid==exon.attrs->Get(i)->attr_id) { |
1251 |
|
|
if (attrval==NULL || |
1252 |
|
|
strcmp(attrval, exon.attrs->Get(i)->attr_val)==0) { |
1253 |
|
|
delcount++; |
1254 |
|
|
exon.attrs->freeItem(i); |
1255 |
|
|
} |
1256 |
|
|
} |
1257 |
|
|
} |
1258 |
|
|
if (delcount>0) exon.attrs->Pack(); |
1259 |
|
|
return delcount; |
1260 |
|
|
} |
1261 |
|
|
|
1262 |
|
|
int GffObj::removeExonAttr(GffExon& exon, int aid, const char* attrval) { |
1263 |
|
|
if (exon.attrs==NULL || aid<0) return 0; |
1264 |
|
|
int delcount=0; //could be more than one |
1265 |
|
|
for (int i=0;i<exon.attrs->Count();i++) { |
1266 |
|
|
if (aid==exon.attrs->Get(i)->attr_id) { |
1267 |
|
|
if (attrval==NULL || |
1268 |
|
|
strcmp(attrval, exon.attrs->Get(i)->attr_val)==0) { |
1269 |
|
|
delcount++; |
1270 |
|
|
exon.attrs->freeItem(i); |
1271 |
|
|
} |
1272 |
|
|
} |
1273 |
|
|
} |
1274 |
|
|
if (delcount>0) exon.attrs->Pack(); |
1275 |
|
|
return delcount; |
1276 |
|
|
} |
1277 |
|
|
|
1278 |
|
|
|
1279 |
|
|
void GffObj::getCDS_ends(uint& cds_start, uint& cds_end) { |
1280 |
|
|
cds_start=0; |
1281 |
|
|
cds_end=0; |
1282 |
|
|
if (CDstart==0 || CDend==0) return; //no CDS info |
1283 |
|
|
int cdsadj=0; |
1284 |
|
|
if (CDphase=='1' || CDphase=='2') { |
1285 |
|
|
cdsadj=CDphase-'0'; |
1286 |
|
|
} |
1287 |
|
|
cds_start=CDstart; |
1288 |
|
|
cds_end=CDend; |
1289 |
|
|
if (strand=='-') cds_end-=cdsadj; |
1290 |
|
|
else cds_start+=cdsadj; |
1291 |
|
|
} |
1292 |
|
|
|
1293 |
|
|
void GffObj::mRNA_CDS_coords(uint& cds_mstart, uint& cds_mend) { |
1294 |
|
|
//sets cds_start and cds_end to the CDS start,end coordinates on the spliced mRNA transcript |
1295 |
|
|
cds_mstart=0; |
1296 |
|
|
cds_mend=0; |
1297 |
|
|
if (CDstart==0 || CDend==0) return; //no CDS info |
1298 |
|
|
//restore normal coordinates, just in case |
1299 |
|
|
unxcoord(); |
1300 |
|
|
int cdsadj=0; |
1301 |
|
|
if (CDphase=='1' || CDphase=='2') { |
1302 |
|
|
cdsadj=CDphase-'0'; |
1303 |
|
|
} |
1304 |
|
|
/* |
1305 |
|
|
uint seqstart=CDstart; |
1306 |
|
|
uint seqend=CDend; |
1307 |
|
|
*/ |
1308 |
|
|
uint seqstart=exons.First()->start; |
1309 |
|
|
uint seqend=exons.Last()->end; |
1310 |
|
|
int s=0; //resulting nucleotide counter |
1311 |
|
|
if (strand=='-') { |
1312 |
|
|
for (int x=exons.Count()-1;x>=0;x--) { |
1313 |
|
|
uint sgstart=exons[x]->start; |
1314 |
|
|
uint sgend=exons[x]->end; |
1315 |
|
|
if (seqend<sgstart || seqstart>sgend) continue; |
1316 |
|
|
if (seqstart>=sgstart && seqstart<=sgend) |
1317 |
|
|
sgstart=seqstart; //seqstart within this segment |
1318 |
|
|
if (seqend>=sgstart && seqend<=sgend) |
1319 |
|
|
sgend=seqend; //seqend within this segment |
1320 |
|
|
s+=(int)(sgend-sgstart)+1; |
1321 |
|
|
if (CDstart>=sgstart && CDstart<=sgend) { |
1322 |
|
|
//CDstart in this segment |
1323 |
|
|
//and we are getting the whole transcript |
1324 |
|
|
cds_mend=s-(int)(CDstart-sgstart); |
1325 |
|
|
} |
1326 |
|
|
if (CDend>=sgstart && CDend<=sgend) { |
1327 |
|
|
//CDstart in this segment |
1328 |
|
|
//and we are getting the whole transcript |
1329 |
|
|
cds_mstart=s-(int)(CDend-cdsadj-sgstart); |
1330 |
|
|
} |
1331 |
|
|
} //for each exon |
1332 |
|
|
} // - strand |
1333 |
|
|
else { // + strand |
1334 |
|
|
for (int x=0;x<exons.Count();x++) { |
1335 |
|
|
uint sgstart=exons[x]->start; |
1336 |
|
|
uint sgend=exons[x]->end; |
1337 |
|
|
if (seqend<sgstart || seqstart>sgend) continue; |
1338 |
|
|
if (seqstart>=sgstart && seqstart<=sgend) |
1339 |
|
|
sgstart=seqstart; //seqstart within this segment |
1340 |
|
|
if (seqend>=sgstart && seqend<=sgend) |
1341 |
|
|
sgend=seqend; //seqend within this segment |
1342 |
|
|
s+=(int)(sgend-sgstart)+1; |
1343 |
|
|
/* for (uint i=sgstart;i<=sgend;i++) { |
1344 |
|
|
spliced[s]=gsubseq[i-gstart]; |
1345 |
|
|
s++; |
1346 |
|
|
}//for each nt |
1347 |
|
|
*/ |
1348 |
|
|
if (CDstart>=sgstart && CDstart<=sgend) { |
1349 |
|
|
//CDstart in this segment |
1350 |
|
|
cds_mstart=s-(int)(sgend-CDstart-cdsadj); |
1351 |
|
|
} |
1352 |
|
|
if (CDend>=sgstart && CDend<=sgend) { |
1353 |
|
|
//CDend in this segment |
1354 |
|
|
cds_mend=s-(int)(sgend-CDend); |
1355 |
|
|
} |
1356 |
|
|
} //for each exon |
1357 |
|
|
} // + strand |
1358 |
|
|
//spliced[s]=0; |
1359 |
|
|
//if (rlen!=NULL) *rlen=s; |
1360 |
|
|
//return spliced; |
1361 |
|
|
} |
1362 |
|
|
|
1363 |
|
|
char* GffObj::getUnspliced(GFaSeqGet* faseq, int* rlen, GList<GSeg>* seglst) |
1364 |
|
|
{ |
1365 |
|
|
if (faseq==NULL) { GMessage("Warning: getUnspliced(NULL,.. ) called!\n"); |
1366 |
|
|
return NULL; |
1367 |
|
|
} |
1368 |
|
|
//restore normal coordinates: |
1369 |
|
|
unxcoord(); |
1370 |
|
|
if (exons.Count()==0) return NULL; |
1371 |
|
|
int fspan=end-start+1; |
1372 |
|
|
const char* gsubseq=faseq->subseq(start, fspan); |
1373 |
|
|
if (gsubseq==NULL) { |
1374 |
|
|
GError("Error getting subseq for %s (%d..%d)!\n", gffID, start, end); |
1375 |
|
|
} |
1376 |
|
|
char* unspliced=NULL; |
1377 |
|
|
|
1378 |
|
|
int seqstart=exons.First()->start; |
1379 |
|
|
int seqend=exons.Last()->end; |
1380 |
|
|
|
1381 |
|
|
int unsplicedlen = 0; |
1382 |
|
|
|
1383 |
|
|
unsplicedlen += seqend - seqstart + 1; |
1384 |
|
|
|
1385 |
|
|
GMALLOC(unspliced, unsplicedlen+1); //allocate more here |
1386 |
|
|
//uint seqstart, seqend; |
1387 |
|
|
|
1388 |
|
|
int s = 0; //resulting nucleotide counter |
1389 |
|
|
if (strand=='-') |
1390 |
|
|
{ |
1391 |
|
|
if (seglst!=NULL) |
1392 |
|
|
seglst->Add(new GSeg(s+1,s+1+seqend-seqstart)); |
1393 |
|
|
for (int i=seqend;i>=seqstart;i--) |
1394 |
|
|
{ |
1395 |
|
|
unspliced[s] = ntComplement(gsubseq[i-start]); |
1396 |
|
|
s++; |
1397 |
|
|
}//for each nt |
1398 |
|
|
} // - strand |
1399 |
|
|
else |
1400 |
|
|
{ // + strand |
1401 |
|
|
if (seglst!=NULL) |
1402 |
|
|
seglst->Add(new GSeg(s+1,s+1+seqend-seqstart)); |
1403 |
|
|
for (int i=seqstart;i<=seqend;i++) |
1404 |
|
|
{ |
1405 |
|
|
unspliced[s]=gsubseq[i-start]; |
1406 |
|
|
s++; |
1407 |
|
|
}//for each nt |
1408 |
|
|
} // + strand |
1409 |
|
|
//assert(s <= unsplicedlen); |
1410 |
|
|
unspliced[s]=0; |
1411 |
|
|
if (rlen!=NULL) *rlen=s; |
1412 |
|
|
return unspliced; |
1413 |
|
|
} |
1414 |
|
|
|
1415 |
|
|
char* GffObj::getSpliced(GFaSeqGet* faseq, bool CDSonly, int* rlen, uint* cds_start, uint* cds_end, |
1416 |
|
|
GList<GSeg>* seglst) { |
1417 |
|
|
if (CDSonly && CDstart==0) return NULL; |
1418 |
|
|
if (faseq==NULL) { GMessage("Warning: getSpliced(NULL,.. ) called!\n"); |
1419 |
|
|
return NULL; |
1420 |
|
|
} |
1421 |
|
|
//restore normal coordinates: |
1422 |
|
|
unxcoord(); |
1423 |
|
|
if (exons.Count()==0) return NULL; |
1424 |
|
|
int fspan=end-start+1; |
1425 |
|
|
const char* gsubseq=faseq->subseq(start, fspan); |
1426 |
|
|
if (gsubseq==NULL) { |
1427 |
|
|
GError("Error getting subseq for %s (%d..%d)!\n", gffID, start, end); |
1428 |
|
|
} |
1429 |
|
|
if (fspan<(int)(end-start+1)) { //special case: stop coordinate was extended past the gseq length, must adjust |
1430 |
|
|
int endadj=end-start+1-fspan; |
1431 |
|
|
uint prevend=end; |
1432 |
|
|
end-=endadj; |
1433 |
|
|
if (CDend>end) CDend=end; |
1434 |
|
|
if (exons.Last()->end>end) { |
1435 |
|
|
exons.Last()->end=end; //this could get us into trouble if exon start is also > end |
1436 |
|
|
if (exons.Last()->start>exons.Last()->end) { |
1437 |
|
|
GError("GffObj::getSpliced() error: improper genomic coordinate %d on %s for %s\n", |
1438 |
|
|
prevend,getGSeqName(), getID()); |
1439 |
|
|
} |
1440 |
|
|
covlen-=endadj; |
1441 |
|
|
} |
1442 |
|
|
} |
1443 |
|
|
char* spliced=NULL; |
1444 |
|
|
GMALLOC(spliced, covlen+1); //allocate more here |
1445 |
|
|
uint seqstart, seqend; |
1446 |
|
|
int cdsadj=0; |
1447 |
|
|
if (CDphase=='1' || CDphase=='2') { |
1448 |
|
|
cdsadj=CDphase-'0'; |
1449 |
|
|
} |
1450 |
|
|
if (CDSonly) { |
1451 |
|
|
seqstart=CDstart; |
1452 |
|
|
seqend=CDend; |
1453 |
|
|
if (strand=='-') seqend-=cdsadj; |
1454 |
|
|
else seqstart+=cdsadj; |
1455 |
|
|
} |
1456 |
|
|
else { |
1457 |
|
|
seqstart=exons.First()->start; |
1458 |
|
|
seqend=exons.Last()->end; |
1459 |
|
|
} |
1460 |
|
|
int s=0; //resulting nucleotide counter |
1461 |
|
|
if (strand=='-') { |
1462 |
|
|
for (int x=exons.Count()-1;x>=0;x--) { |
1463 |
|
|
uint sgstart=exons[x]->start; |
1464 |
|
|
uint sgend=exons[x]->end; |
1465 |
|
|
if (seqend<sgstart || seqstart>sgend) continue; |
1466 |
|
|
if (seqstart>=sgstart && seqstart<=sgend) |
1467 |
|
|
sgstart=seqstart; //seqstart within this segment |
1468 |
|
|
if (seqend>=sgstart && seqend<=sgend) |
1469 |
|
|
sgend=seqend; //seqend within this segment |
1470 |
|
|
if (seglst!=NULL) |
1471 |
|
|
seglst->Add(new GSeg(s+1,s+1+sgend-sgstart)); |
1472 |
|
|
for (uint i=sgend;i>=sgstart;i--) { |
1473 |
|
|
spliced[s] = ntComplement(gsubseq[i-start]); |
1474 |
|
|
s++; |
1475 |
|
|
}//for each nt |
1476 |
|
|
|
1477 |
|
|
if (!CDSonly && cds_start!=NULL && CDstart>0) { |
1478 |
|
|
if (CDstart>=sgstart && CDstart<=sgend) { |
1479 |
|
|
//CDstart in this segment |
1480 |
|
|
//and we are getting the whole transcript |
1481 |
|
|
*cds_end=s-(CDstart-sgstart); |
1482 |
|
|
} |
1483 |
|
|
if (CDend>=sgstart && CDend<=sgend) { |
1484 |
|
|
//CDstart in this segment |
1485 |
|
|
//and we are getting the whole transcript |
1486 |
|
|
*cds_start=s-(CDend-cdsadj-sgstart); |
1487 |
|
|
} |
1488 |
|
|
}//update local CDS coordinates |
1489 |
|
|
} //for each exon |
1490 |
|
|
} // - strand |
1491 |
|
|
else { // + strand |
1492 |
|
|
for (int x=0;x<exons.Count();x++) { |
1493 |
|
|
uint sgstart=exons[x]->start; |
1494 |
|
|
uint sgend=exons[x]->end; |
1495 |
|
|
if (seqend<sgstart || seqstart>sgend) continue; |
1496 |
|
|
if (seqstart>=sgstart && seqstart<=sgend) |
1497 |
|
|
sgstart=seqstart; //seqstart within this segment |
1498 |
|
|
if (seqend>=sgstart && seqend<=sgend) |
1499 |
|
|
sgend=seqend; //seqend within this segment |
1500 |
|
|
if (seglst!=NULL) |
1501 |
|
|
seglst->Add(new GSeg(s+1,s+1+sgend-sgstart)); |
1502 |
|
|
for (uint i=sgstart;i<=sgend;i++) { |
1503 |
|
|
spliced[s]=gsubseq[i-start]; |
1504 |
|
|
s++; |
1505 |
|
|
}//for each nt |
1506 |
|
|
if (!CDSonly && cds_start!=NULL && CDstart>0) { |
1507 |
|
|
if (CDstart>=sgstart && CDstart<=sgend) { |
1508 |
|
|
//CDstart in this segment |
1509 |
|
|
//and we are getting the whole transcript |
1510 |
|
|
*cds_start=s-(sgend-CDstart-cdsadj); |
1511 |
|
|
} |
1512 |
|
|
if (CDend>=sgstart && CDend<=sgend) { |
1513 |
|
|
//CDstart in this segment |
1514 |
|
|
//and we are getting the whole transcript |
1515 |
|
|
*cds_end=s-(sgend-CDend); |
1516 |
|
|
} |
1517 |
|
|
}//update local CDS coordinates |
1518 |
|
|
} //for each exon |
1519 |
|
|
} // + strand |
1520 |
|
|
spliced[s]=0; |
1521 |
|
|
if (rlen!=NULL) *rlen=s; |
1522 |
|
|
return spliced; |
1523 |
|
|
} |
1524 |
|
|
|
1525 |
|
|
char* GffObj::getSplicedTr(GFaSeqGet* faseq, bool CDSonly, int* rlen) { |
1526 |
|
|
if (CDSonly && CDstart==0) return NULL; |
1527 |
|
|
//restore normal coordinates: |
1528 |
|
|
unxcoord(); |
1529 |
|
|
if (exons.Count()==0) return NULL; |
1530 |
|
|
int fspan=end-start+1; |
1531 |
|
|
const char* gsubseq=faseq->subseq(start, fspan); |
1532 |
|
|
if (gsubseq==NULL) { |
1533 |
|
|
GError("Error getting subseq for %s (%d..%d)!\n", gffID, start, end); |
1534 |
|
|
} |
1535 |
|
|
|
1536 |
|
|
char* translation=NULL; |
1537 |
|
|
GMALLOC(translation, (int)(covlen/3)+1); |
1538 |
|
|
uint seqstart, seqend; |
1539 |
|
|
int cdsadj=0; |
1540 |
|
|
if (CDphase=='1' || CDphase=='2') { |
1541 |
|
|
cdsadj=CDphase-'0'; |
1542 |
|
|
} |
1543 |
|
|
if (CDSonly) { |
1544 |
|
|
seqstart=CDstart; |
1545 |
|
|
seqend=CDend; |
1546 |
|
|
if (strand=='-') seqend-=cdsadj; |
1547 |
|
|
else seqstart+=cdsadj; |
1548 |
|
|
} |
1549 |
|
|
else { |
1550 |
|
|
seqstart=exons.First()->start; |
1551 |
|
|
seqend=exons.Last()->end; |
1552 |
|
|
} |
1553 |
|
|
Codon codon; |
1554 |
|
|
int nt=0; //codon nucleotide counter (0..2) |
1555 |
|
|
int aa=0; //aminoacid count |
1556 |
|
|
if (strand=='-') { |
1557 |
|
|
for (int x=exons.Count()-1;x>=0;x--) { |
1558 |
|
|
uint sgstart=exons[x]->start; |
1559 |
|
|
uint sgend=exons[x]->end; |
1560 |
|
|
if (seqend<sgstart || seqstart>sgend) continue; |
1561 |
|
|
if (seqstart>=sgstart && seqstart<=sgend) |
1562 |
|
|
sgstart=seqstart; //seqstart within this segment |
1563 |
|
|
if (seqend>=sgstart && seqend<=sgend) { |
1564 |
|
|
sgend=seqend; //seqend within this segment |
1565 |
|
|
} |
1566 |
|
|
for (uint i=sgend;i>=sgstart;i--) { |
1567 |
|
|
codon.nuc[nt]=ntComplement(gsubseq[i-start]); |
1568 |
|
|
nt++; |
1569 |
|
|
if (nt==3) { |
1570 |
|
|
nt=0; |
1571 |
|
|
translation[aa]=codon.translate(); |
1572 |
|
|
aa++; |
1573 |
|
|
} |
1574 |
|
|
}//for each nt |
1575 |
|
|
} //for each exon |
1576 |
|
|
} // - strand |
1577 |
|
|
else { // + strand |
1578 |
|
|
for (int x=0;x<exons.Count();x++) { |
1579 |
|
|
uint sgstart=exons[x]->start; |
1580 |
|
|
uint sgend=exons[x]->end; |
1581 |
|
|
if (seqend<sgstart || seqstart>sgend) continue; |
1582 |
|
|
if (seqstart>=sgstart && seqstart<=sgend) |
1583 |
|
|
sgstart=seqstart; //seqstart within this segment |
1584 |
|
|
if (seqend>=sgstart && seqend<=sgend) |
1585 |
|
|
sgend=seqend; //seqend within this segment |
1586 |
|
|
for (uint i=sgstart;i<=sgend;i++) { |
1587 |
|
|
codon.nuc[nt]=gsubseq[i-start]; |
1588 |
|
|
nt++; |
1589 |
|
|
if (nt==3) { |
1590 |
|
|
nt=0; |
1591 |
|
|
translation[aa]=codon.translate(); |
1592 |
|
|
aa++; |
1593 |
|
|
} |
1594 |
|
|
}//for each nt |
1595 |
|
|
} //for each exon |
1596 |
|
|
} // + strand |
1597 |
|
|
translation[aa]=0; |
1598 |
|
|
if (rlen!=NULL) *rlen=aa; |
1599 |
|
|
return translation; |
1600 |
|
|
} |
1601 |
|
|
|
1602 |
|
|
void GffObj::printSummary(FILE* fout) { |
1603 |
|
|
if (fout==NULL) fout=stdout; |
1604 |
|
|
fprintf(fout, "%s\t%c\t%d\t%d\t%4.2f\t%4.1f\n", gffID, |
1605 |
|
|
strand, start, end, gscore, (float)qcov/10.0); |
1606 |
|
|
} |
1607 |
|
|
|
1608 |
|
|
void GffObj::printGxfLine(FILE* fout, const char* tlabel, const char* gseqname, bool iscds, |
1609 |
|
|
uint segstart, uint segend, int exidx, char phase, bool gff3) { |
1610 |
|
|
static char scorestr[14]; |
1611 |
|
|
strcpy(scorestr,"."); |
1612 |
|
|
GffAttrs* xattrs=NULL; |
1613 |
|
|
if (exidx>=0) { |
1614 |
|
|
if (exons[exidx]->score) sprintf(scorestr,"%.2f", exons[exidx]->score); |
1615 |
|
|
xattrs=exons[exidx]->attrs; |
1616 |
|
|
} |
1617 |
|
|
if (phase==0 || !iscds) phase='.'; |
1618 |
|
|
const char* ftype=iscds ? "CDS" : getSubfName(); |
1619 |
|
|
if (gff3) { |
1620 |
|
|
fprintf(fout, |
1621 |
|
|
"%s\t%s\t%s\t%d\t%d\t%s\t%c\t%c\tParent=%s", |
1622 |
|
|
gseqname, tlabel, ftype, segstart, segend, scorestr, strand, |
1623 |
|
|
phase, gffID); |
1624 |
|
|
if (xattrs!=NULL) { |
1625 |
|
|
for (int i=0;i<xattrs->Count();i++) |
1626 |
|
|
fprintf(fout, ";%s=%s",names->attrs.getName(xattrs->Get(i)->attr_id), |
1627 |
|
|
xattrs->Get(i)->attr_val); |
1628 |
|
|
} |
1629 |
|
|
fprintf(fout, "\n"); |
1630 |
|
|
} //GFF |
1631 |
|
|
else {//for GTF -- we print only transcripts |
1632 |
|
|
//if (isValidTranscript()) |
1633 |
|
|
fprintf(fout, "%s\t%s\t%s\t%d\t%d\t%s\t%c\t%c\ttranscript_id \"%s\";", |
1634 |
|
|
gseqname, tlabel, ftype, segstart, segend, scorestr, strand, phase, gffID); |
1635 |
|
|
//char* geneid=(geneID!=NULL)? geneID : gffID; |
1636 |
|
|
if (geneID) |
1637 |
|
|
fprintf(fout," gene_id \"%s\";",geneID); |
1638 |
|
|
if (gene_name!=NULL) { |
1639 |
|
|
//fprintf(fout, " gene_name "); |
1640 |
|
|
//if (gene_name[0]=='"') fprintf (fout, "%s;",gene_name); |
1641 |
|
|
// else fprintf(fout, "\"%s\";",gene_name); |
1642 |
|
|
fprintf(fout," gene_name \"%s\";",gene_name); |
1643 |
|
|
} |
1644 |
|
|
if (xattrs!=NULL) { |
1645 |
|
|
for (int i=0;i<xattrs->Count();i++) { |
1646 |
|
|
if (xattrs->Get(i)->attr_val==NULL) continue; |
1647 |
|
|
const char* attrname=names->attrs.getName(xattrs->Get(i)->attr_id); |
1648 |
|
|
fprintf(fout, " %s ",attrname); |
1649 |
|
|
if (xattrs->Get(i)->attr_val[0]=='"') |
1650 |
|
|
fprintf(fout, "%s;",xattrs->Get(i)->attr_val); |
1651 |
|
|
else fprintf(fout, "\"%s\";",xattrs->Get(i)->attr_val); |
1652 |
|
|
} |
1653 |
|
|
} |
1654 |
|
|
//for GTF, also append the GffObj attributes to each exon line |
1655 |
|
|
if ((xattrs=this->attrs)!=NULL) { |
1656 |
|
|
for (int i=0;i<xattrs->Count();i++) { |
1657 |
|
|
if (xattrs->Get(i)->attr_val==NULL) continue; |
1658 |
|
|
const char* attrname=names->attrs.getName(xattrs->Get(i)->attr_id); |
1659 |
|
|
fprintf(fout, " %s ",attrname); |
1660 |
|
|
if (xattrs->Get(i)->attr_val[0]=='"') |
1661 |
|
|
fprintf(fout, "%s;",xattrs->Get(i)->attr_val); |
1662 |
|
|
else fprintf(fout, "\"%s\";",xattrs->Get(i)->attr_val); |
1663 |
|
|
} |
1664 |
|
|
} |
1665 |
|
|
fprintf(fout, "\n"); |
1666 |
|
|
}//GTF |
1667 |
|
|
} |
1668 |
|
|
|
1669 |
|
|
void GffObj::printGxf(FILE* fout, GffPrintMode gffp, |
1670 |
|
|
const char* tlabel, const char* gfparent) { |
1671 |
|
|
static char tmpstr[255]; |
1672 |
|
|
if (tlabel==NULL) { |
1673 |
|
|
tlabel=track_id>=0 ? names->tracks.Get(track_id)->name : |
1674 |
|
|
(char*)"gffobj" ; |
1675 |
|
|
} |
1676 |
|
|
unxcoord(); |
1677 |
|
|
//if (exons.Count()==0) return; |
1678 |
|
|
const char* gseqname=names->gseqs.Get(gseq_id)->name; |
1679 |
|
|
bool gff3 = (gffp>=pgffAny); |
1680 |
|
|
bool showCDS = (gffp==pgtfAny || gffp==pgtfCDS || gffp==pgffCDS || gffp==pgffAny || gffp==pgffBoth); |
1681 |
|
|
bool showExon = (gffp<=pgtfExon || gffp==pgffAny || gffp==pgffExon || gffp==pgffBoth); |
1682 |
|
|
if (gff3) { |
1683 |
|
|
//print GFF3 mRNA line: |
1684 |
|
|
if (gscore>0.0) sprintf(tmpstr,"%.2f", gscore); |
1685 |
|
|
else strcpy(tmpstr,"."); |
1686 |
|
|
uint pstart, pend; |
1687 |
|
|
if (gffp==pgffCDS) { |
1688 |
|
|
pstart=CDstart; |
1689 |
|
|
pend=CDend; |
1690 |
|
|
} |
1691 |
|
|
else { pstart=start;pend=end; } |
1692 |
|
|
//const char* ftype=isTranscript() ? "mRNA" : getFeatureName(); |
1693 |
|
|
const char* ftype=getFeatureName(); |
1694 |
|
|
fprintf(fout, |
1695 |
|
|
"%s\t%s\t%s\t%d\t%d\t%s\t%c\t.\tID=%s", |
1696 |
|
|
gseqname, tlabel, ftype, pstart, pend, tmpstr, strand, gffID); |
1697 |
|
|
if (CDstart>0 && !showCDS && !isCDS) fprintf(fout,";CDS=%d-%d",CDstart,CDend); |
1698 |
|
|
if (gfparent!=NULL) { |
1699 |
|
|
//parent override |
1700 |
|
|
fprintf(fout, ";Parent=%s",gfparent); |
1701 |
|
|
} |
1702 |
|
|
else { |
1703 |
|
|
if (parent!=NULL && !parent->isDiscarded()) |
1704 |
|
|
fprintf(fout, ";Parent=%s",parent->getID()); |
1705 |
|
|
} |
1706 |
|
|
if (geneID!=NULL) |
1707 |
|
|
fprintf(fout, ";geneID=%s",geneID); |
1708 |
|
|
if (gene_name!=NULL) |
1709 |
|
|
fprintf(fout, ";gene_name=%s",gene_name); |
1710 |
|
|
if (attrs!=NULL) { |
1711 |
|
|
for (int i=0;i<attrs->Count();i++) { |
1712 |
|
|
const char* attrname=names->attrs.getName(attrs->Get(i)->attr_id); |
1713 |
|
|
fprintf(fout,";%s=%s", attrname, |
1714 |
|
|
attrs->Get(i)->attr_val); |
1715 |
|
|
} |
1716 |
|
|
} |
1717 |
|
|
fprintf(fout,"\n"); |
1718 |
|
|
}// gff3 mRNA line |
1719 |
|
|
if (showExon) { |
1720 |
|
|
//print exons |
1721 |
|
|
if (isCDS && exons.Count()>0 && |
1722 |
|
|
((strand=='-' && exons.Last()->phase<'0') || (strand=='+' && exons.Last()->phase<'0'))) |
1723 |
|
|
updateExonPhase(); |
1724 |
|
|
|
1725 |
|
|
for (int i=0;i<exons.Count();i++) { |
1726 |
|
|
printGxfLine(fout, tlabel, gseqname, isCDS, exons[i]->start, exons[i]->end, i, exons[i]->phase, gff3); |
1727 |
|
|
} |
1728 |
|
|
}//printing exons |
1729 |
|
|
if (showCDS && !isCDS && CDstart>0) { |
1730 |
|
|
GArray<GffCDSeg> cds(true,true); |
1731 |
|
|
getCDSegs(cds); |
1732 |
|
|
for (int i=0;i<cds.Count();i++) { |
1733 |
|
|
printGxfLine(fout, tlabel, gseqname, true, cds[i].start, cds[i].end, -1, cds[i].phase, gff3); |
1734 |
|
|
} |
1735 |
|
|
} //showCDS |
1736 |
|
|
} |
1737 |
|
|
|
1738 |
|
|
void GffObj::updateExonPhase() { |
1739 |
|
|
if (!isCDS) return; |
1740 |
|
|
int cdsacc=0; |
1741 |
|
|
if (CDphase=='1' || CDphase=='2') { |
1742 |
|
|
cdsacc+= 3-(CDphase-'0'); |
1743 |
|
|
} |
1744 |
|
|
if (strand=='-') { //reverse strand |
1745 |
|
|
for (int i=exons.Count()-1;i>=0;i--) { |
1746 |
|
|
exons[i]->phase='0'+ (3-cdsacc%3)%3; |
1747 |
|
|
cdsacc+=exons[i]->end-exons[i]->start+1; |
1748 |
|
|
} |
1749 |
|
|
} |
1750 |
|
|
else { //forward strand |
1751 |
|
|
for (int i=0;i<exons.Count();i++) { |
1752 |
|
|
exons[i]->phase='0'+ (3-cdsacc%3)%3; |
1753 |
|
|
cdsacc+=exons[i]->end-exons[i]->start+1; |
1754 |
|
|
} |
1755 |
|
|
} |
1756 |
|
|
} |
1757 |
|
|
|
1758 |
|
|
|
1759 |
|
|
void GffObj::getCDSegs(GArray<GffCDSeg>& cds) { |
1760 |
|
|
GffCDSeg cdseg; |
1761 |
|
|
int cdsacc=0; |
1762 |
|
|
if (CDphase=='1' || CDphase=='2') { |
1763 |
|
|
cdsacc+= 3-(CDphase-'0'); |
1764 |
|
|
} |
1765 |
|
|
if (strand=='-') { |
1766 |
|
|
for (int x=exons.Count()-1;x>=0;x--) { |
1767 |
|
|
uint sgstart=exons[x]->start; |
1768 |
|
|
uint sgend=exons[x]->end; |
1769 |
|
|
if (CDend<sgstart || CDstart>sgend) continue; |
1770 |
|
|
if (CDstart>=sgstart && CDstart<=sgend) |
1771 |
|
|
sgstart=CDstart; //cdstart within this segment |
1772 |
|
|
if (CDend>=sgstart && CDend<=sgend) |
1773 |
|
|
sgend=CDend; //cdend within this segment |
1774 |
|
|
cdseg.start=sgstart; |
1775 |
|
|
cdseg.end=sgend; |
1776 |
|
|
cdseg.exonidx=x; |
1777 |
|
|
//cdseg.phase='0'+(cdsacc>0 ? (3-cdsacc%3)%3 : 0); |
1778 |
|
|
cdseg.phase='0'+ (3-cdsacc%3)%3; |
1779 |
|
|
cdsacc+=sgend-sgstart+1; |
1780 |
|
|
cds.Add(cdseg); |
1781 |
|
|
} //for each exon |
1782 |
|
|
} // - strand |
1783 |
|
|
else { // + strand |
1784 |
|
|
for (int x=0;x<exons.Count();x++) { |
1785 |
|
|
uint sgstart=exons[x]->start; |
1786 |
|
|
uint sgend=exons[x]->end; |
1787 |
|
|
if (CDend<sgstart || CDstart>sgend) continue; |
1788 |
|
|
if (CDstart>=sgstart && CDstart<=sgend) |
1789 |
|
|
sgstart=CDstart; //seqstart within this segment |
1790 |
|
|
if (CDend>=sgstart && CDend<=sgend) |
1791 |
|
|
sgend=CDend; //seqend within this segment |
1792 |
|
|
cdseg.start=sgstart; |
1793 |
|
|
cdseg.end=sgend; |
1794 |
|
|
cdseg.exonidx=x; |
1795 |
|
|
//cdseg.phase='0'+(cdsacc>0 ? (3-cdsacc%3)%3 : 0); |
1796 |
|
|
cdseg.phase='0' + (3-cdsacc%3)%3 ; |
1797 |
|
|
cdsacc+=sgend-sgstart+1; |
1798 |
|
|
cds.Add(cdseg); |
1799 |
|
|
} //for each exon |
1800 |
|
|
} // + strand |
1801 |
|
|
} |