ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/gclib/gff.h
Revision: 227
Committed: Mon Mar 26 21:28:07 2012 UTC (7 years, 4 months ago) by gpertea
File size: 35271 byte(s)
Log Message:
Line User Rev File contents
1 gpertea 2 #ifndef GFF_H
2     #define GFF_H
3    
4     #include "GBase.h"
5     #include "gdna.h"
6     #include "codons.h"
7     #include "GFaSeqGet.h"
8     #include "GList.hh"
9     #include "GHash.hh"
10    
11     /*
12     const byte exMskMajSpliceL = 0x01;
13     const byte exMskMajSpliceR = 0x02;
14     const byte exMskMinSpliceL = 0x04;
15     const byte exMskMinSpliceR = 0x08;
16     const byte exMskTag = 0x80;
17     */
18 gpertea 16
19     //reserved Gffnames::feats entries -- basic feature types
20     extern const int gff_fid_mRNA; // "mRNA" feature name
21     extern const int gff_fid_transcript; // *RNA, *transcript feature name
22 gpertea 2 extern const int gff_fid_exon;
23 gpertea 16 extern const int gff_fid_CDS; //never really used, except for display only
24     //use gff_fid_exon instead
25 gpertea 2 extern const uint GFF_MAX_LOCUS;
26     extern const uint GFF_MAX_EXON;
27     extern const uint GFF_MAX_INTRON;
28    
29 gpertea 16 extern const uint gfo_flag_CHILDREN_PROMOTED;
30     extern const uint gfo_flag_HAS_ERRORS;
31     extern const uint gfo_flag_IS_GENE;
32 gpertea 153 extern const uint gfo_flag_HAS_GFF_ID; //found a GFF3 formatted main feature with its own ID
33 gpertea 16 extern const uint gfo_flag_BY_EXON; //created by subfeature (exon) directly
34     //(GTF2 and some chado gff3 dumps with exons given before their mRNA)
35     extern const uint gfo_flag_IS_TRANSCRIPT; //recognized as '*RNA' or '*transcript'
36     extern const uint gfo_flag_DISCARDED; //should not be printed under the "transcriptsOnly" directive
37     extern const uint gfo_flag_LST_KEEP; //GffObj from GffReader::gflst is to be kept (not deallocated)
38     //when GffReader is destroyed
39     extern const uint gfo_flag_LEVEL_MSK; //hierarchical level: 0 = no parent
40     extern const byte gfo_flagShift_LEVEL;
41    
42     extern bool gff_show_warnings;
43    
44 gpertea 2 #define GFF_LINELEN 2048
45     #define ERR_NULL_GFNAMES "Error: GffObj::%s requires a non-null GffNames* names!\n"
46    
47    
48     enum GffExonType {
49 gpertea 16 exgffNone=0, //not a recognizable exon or CDS segment
50     exgffStart, //from "start_codon" feature (within CDS)
51     exgffStop, //from "stop_codon" feature (may be outside CDS)
52 gpertea 2 exgffCDS, //from "CDS" feature
53     exgffUTR, //from "UTR" feature
54 gpertea 16 exgffCDSUTR, //from a merge of UTR and CDS feature
55 gpertea 2 exgffExon, //from "exon" feature
56     };
57    
58     class GffReader;
59    
60     class GffLine {
61 gpertea 16 char* _parents; //stores a copy of the Parent attribute value,
62     //with commas replaced by \0
63     int _parents_len;
64 gpertea 2 public:
65 gpertea 16 char* dupline; //duplicate of original line
66     char* line; //this will have tabs replaced by \0
67     int llen;
68 gpertea 2 char* gseqname;
69     char* track;
70     char* ftype; //feature name: mRNA/gene/exon/CDS
71 gpertea 16 char* info; //the last, attributes' field, unparsed
72 gpertea 2 uint fstart;
73     uint fend;
74     uint qstart; //overlap coords on query, if available
75     uint qend;
76     uint qlen; //query len, if given
77     double score;
78     char strand;
79     bool skip;
80 gpertea 16 bool is_gff3; //if the line appears to be in GFF3 format
81 gpertea 2 bool is_cds; //"cds" and "stop_codon" features
82     bool is_exon; //"exon" and "utr" features
83     char exontype; // gffExonType
84 gpertea 16 bool is_transcript; //if current feature is *RNA or *transcript
85     bool is_gene; //if current feature is *gene
86 gpertea 2 char phase; // '.' , '0', '1' or '2'
87 gpertea 16 // -- allocated strings:
88     char* gene_name; //value of gene_name attribute (GTF) if present or Name attribute of a gene feature (GFF3)
89     char* gene_id; //value of gene_id attribute (GTF) if present or ID attribute of a gene feature (GFF3)
90 gpertea 2 //
91 gpertea 16 char** parents; //for GTF only parents[0] is used
92     int num_parents;
93     char* ID; // if a ID=.. attribute was parsed, or a GTF with 'transcript' line (transcript_id)
94 gpertea 2 GffLine(GffReader* reader, const char* l); //parse the line accordingly
95 gpertea 16 void discardParent() {
96     GFREE(_parents);
97     _parents_len=0;
98     num_parents=0;
99     parents=NULL;
100     }
101 gpertea 150 char* extractAttr(const char* pre, bool caseStrict=false, bool enforce_GTF2=false);
102 gpertea 16 GffLine(GffLine* l) { //a copy constructor
103     memcpy((void*)this, (void*)l, sizeof(GffLine));
104     line=NULL;
105     GMALLOC(line, llen+1);
106     memcpy(line, l->line, llen+1);
107     GMALLOC(dupline, llen+1);
108     memcpy(dupline, l->dupline, llen+1);
109     //--offsets within line[]
110     gseqname=line+(l->gseqname-l->line);
111     track=line+(l->track-l->line);
112     ftype=line+(l->ftype-l->line);
113     info=line+(l->info-l->line);
114     //Parent=Gstrdup(l->Parent);
115     if (l->_parents_len>0) {
116     _parents_len=l->_parents_len;
117     GMALLOC(_parents, _parents_len);
118     memcpy(_parents, l->_parents, _parents_len);
119     num_parents=l->num_parents;
120     for (int i=0;i<num_parents;i++) {
121     parents[i]=_parents+(l->parents[i] - l->_parents);
122     }
123     }
124     //-- allocated string copies:
125     ID=Gstrdup(l->ID);
126     if (l->gene_name!=NULL)
127     gene_name=Gstrdup(l->gene_name);
128     if (l->gene_id!=NULL)
129     gene_id=Gstrdup(l->gene_id);
130     }
131 gpertea 2 GffLine() {
132     line=NULL;
133 gpertea 16 dupline=NULL;
134 gpertea 2 gseqname=NULL;
135     track=NULL;
136     ftype=NULL;
137     fstart=0;
138     fend=0;
139     info=NULL;
140 gpertea 16 _parents=NULL;
141     _parents_len=0;
142     parents=NULL;
143     num_parents=0;
144 gpertea 2 ID=NULL;
145 gpertea 16 gene_name=NULL;
146     gene_id=NULL;
147 gpertea 2 skip=true;
148     qstart=0;
149     qend=0;
150     qlen=0;
151     exontype=0;
152     is_cds=false;
153 gpertea 16 is_gff3=false;
154     is_transcript=false;
155     is_gene=false;
156 gpertea 2 is_exon=false;
157     }
158     ~GffLine() {
159 gpertea 16 GFREE(dupline);
160 gpertea 2 GFREE(line);
161 gpertea 16 GFREE(_parents);
162     GFREE(parents);
163 gpertea 2 GFREE(ID);
164 gpertea 16 GFREE(gene_name);
165     GFREE(gene_id);
166 gpertea 2 }
167     };
168    
169     class GffAttr {
170     public:
171     int attr_id;
172     char* attr_val;
173     GffAttr(int an_id, const char* av=NULL) {
174     attr_id=an_id;
175     attr_val=NULL;
176 gpertea 16 setValue(av);
177 gpertea 2 }
178     ~GffAttr() {
179     GFREE(attr_val);
180     }
181 gpertea 16 void setValue(const char* av) {
182     if (attr_val!=NULL) {
183     GFREE(attr_val);
184     }
185     if (av==NULL || av[0]==0) return;
186     //trim spaces
187     const char* vstart=av;
188     while (*vstart==' ') av++;
189     const char* vend=vstart;
190     bool keep_dq=false;
191     while (vend[1]!=0) {
192     if (*vend==' ' && vend[1]!=' ') keep_dq=true;
193     else if (*vend==';') keep_dq=true;
194     vend++;
195     }
196     //remove spaces at the end:
197     while (*vend==' ' && vend!=vstart) vend--;
198     //practical clean-up: if it doesn't have any internal spaces just strip those useless double quotes
199     if (!keep_dq && *vstart=='"' && *vend=='"') {
200     vend--;
201     vstart++;
202     }
203     attr_val=Gstrdup(vstart, vend);
204     }
205 gpertea 2 bool operator==(GffAttr& d){
206     return (this==&d);
207     }
208     bool operator>(GffAttr& d){
209     return (this>&d);
210     }
211     bool operator<(GffAttr& d){
212     return (this<&d);
213     }
214    
215     };
216    
217     class GffNameList;
218     class GffNames;
219    
220     class GffNameInfo {
221     friend class GffNameList;
222     protected:
223     int idx;
224     public:
225     char* name;
226     GffNameInfo() { name=NULL; idx=-1; }
227     GffNameInfo(const char* n) {
228     name=Gstrdup(n);
229     }
230    
231     ~GffNameInfo() {
232 gpertea 16 GFREE(name);
233     }
234 gpertea 2
235     bool operator==(GffNameInfo& d){
236     return (strcmp(this->name, d.name)==0);
237     }
238     bool operator<(GffNameInfo& d){
239     return (strcmp(this->name, d.name)<0);
240     }
241     };
242    
243     class GffNameList:public GList<GffNameInfo> {
244     friend class GffNameInfo;
245     friend class GffNames;
246     protected:
247     GHash<GffNameInfo> byName;//hash with shared keys
248     int idlast; //fList index of last added/reused name
249     void addStatic(const char* tname) {// fast add
250     GffNameInfo* f=new GffNameInfo(tname);
251     idlast=this->Add(f);
252     f->idx=idlast;
253     byName.shkAdd(f->name,f);
254     }
255     public:
256 gpertea 227 GffNameList(int init_capacity=6):GList<GffNameInfo>(init_capacity, false,true,true), byName(false) {
257 gpertea 2 idlast=-1;
258 gpertea 227 setCapacity(init_capacity);
259 gpertea 2 }
260     char* lastNameUsed() { return idlast<0 ? NULL : Get(idlast)->name; }
261     int lastNameId() { return idlast; }
262     char* getName(int nid) { //retrieve name by its ID
263     if (nid<0 || nid>=fCount)
264     GError("GffNameList Error: invalid index (%d)\n",nid);
265     return fList[nid]->name;
266     }
267    
268     int addName(const char* tname) {//returns or create an id for the given name
269     //check idlast first, chances are it's the same feature name checked
270 gpertea 227 /*if (idlast>=0 && strcmp(fList[idlast]->name,tname)==0)
271     return idlast;*/
272 gpertea 2 GffNameInfo* f=byName.Find(tname);
273     int fidx=-1;
274     if (f!=NULL) fidx=f->idx;
275     else {//add new entry
276     f=new GffNameInfo(tname);
277     fidx=this->Add(f);
278     f->idx=fidx;
279     byName.shkAdd(f->name,f);
280     }
281     idlast=fidx;
282     return fidx;
283     }
284 gpertea 16
285     int addNewName(const char* tname) {
286     GffNameInfo* f=new GffNameInfo(tname);
287     int fidx=this->Add(f);
288     f->idx=fidx;
289     byName.shkAdd(f->name,f);
290     return fidx;
291     }
292    
293 gpertea 2 int getId(const char* tname) { //only returns a name id# if found
294     GffNameInfo* f=byName.Find(tname);
295     if (f==NULL) return -1;
296     return f->idx;
297     }
298 gpertea 16 int removeName() {
299 gpertea 2 GError("Error: removing names from GffNameList not allowed!\n");
300     return -1;
301     }
302     };
303    
304     class GffNames {
305     public:
306     int numrefs;
307     GffNameList tracks;
308     GffNameList gseqs;
309     GffNameList attrs;
310 gpertea 16 GffNameList feats; //feature names: 'mRNA', 'exon', 'CDS' etc.
311 gpertea 2 GffNames():tracks(),gseqs(),attrs(), feats() {
312     numrefs=0;
313     //the order below is critical!
314     //has to match: gff_fid_mRNA, gff_fid_exon, gff_fid_CDS
315     feats.addStatic("mRNA");//index 0=gff_fid_mRNA
316 gpertea 16 feats.addStatic("transcript");//index 1=gff_fid_transcript
317 gpertea 2 feats.addStatic("exon");//index 1=gff_fid_exon
318     feats.addStatic("CDS"); //index 2=gff_fid_CDS
319     }
320     };
321    
322     void gffnames_ref(GffNames* &n);
323     void gffnames_unref(GffNames* &n);
324    
325     enum GffPrintMode {
326     pgtfAny, //print record as read
327     pgtfExon,
328     pgtfCDS,
329     pgffAny, //print record as read
330     pgffExon,
331     pgffCDS,
332     pgffBoth,
333     };
334    
335    
336     class GffAttrs:public GList<GffAttr> {
337     public:
338     GffAttrs():GList<GffAttr>(false,true,false) { }
339 gpertea 16 void add_or_update(GffNames* names, const char* attrname, const char* val) {
340     int aid=names->attrs.getId(attrname);
341     if (aid>=0) {
342     //attribute found in the dictionary
343     for (int i=0;i<Count();i++) {
344     //do we have it?
345     if (aid==Get(i)->attr_id) {
346     //update the value
347     Get(i)->setValue(val);
348     return;
349     }
350     }
351     }
352     else {
353     aid=names->attrs.addNewName(attrname);
354     }
355     this->Add(new GffAttr(aid, val));
356     }
357    
358 gpertea 2 char* getAttr(GffNames* names, const char* attrname) {
359     int aid=names->attrs.getId(attrname);
360     if (aid>=0)
361     for (int i=0;i<Count();i++)
362     if (aid==Get(i)->attr_id) return Get(i)->attr_val;
363     return NULL;
364     }
365 gpertea 16 char* getAttr(int aid) {
366     if (aid>=0)
367     for (int i=0;i<Count();i++)
368     if (aid==Get(i)->attr_id) return Get(i)->attr_val;
369     return NULL;
370     }
371 gpertea 2 };
372    
373    
374     class GffExon : public GSeg {
375     public:
376 gpertea 16 void* uptr; //for later extensions
377 gpertea 2 GffAttrs* attrs; //other attributes kept for this exon
378     double score; // gff score column
379     char phase; //GFF phase column - for CDS segments only
380     // '.' = undefined (UTR), '0','1','2' for CDS exons
381     char exontype; // 1="exon" 2="cds" 3="utr" 4="stop_codon"
382     int qstart; // for mRNA/protein exon mappings: coordinates on query
383     int qend;
384     GffExon(int s=0, int e=0, double sc=0, char fr=0, int qs=0, int qe=0, char et=0) {
385     uptr=NULL;
386     attrs=NULL;
387     if (s<e) {
388     start=s;
389     end=e;
390     }
391     else {
392     start=e;
393     end=s;
394     }
395     if (qs<qe) {
396     qstart=qs;
397     qend=qe;
398     } else {
399     qstart=qe;
400     qend=qs;
401     }
402     score=sc;
403     phase=fr;
404     exontype=et;
405     } //constructor
406    
407     char* getAttr(GffNames* names, const char* atrname) {
408     if (attrs==NULL || names==NULL || atrname==NULL) return NULL;
409     return attrs->getAttr(names, atrname);
410     }
411    
412 gpertea 16 char* getAttr(int aid) {
413     if (attrs==NULL) return NULL;
414     return attrs->getAttr(aid);
415     }
416    
417 gpertea 2 ~GffExon() { //destructor
418     if (attrs!=NULL) delete attrs;
419     }
420     };
421    
422    
423     class GffCDSeg:public GSeg {
424     public:
425     char phase;
426     int exonidx;
427     };
428     //one GFF mRNA object -- e.g. a mRNA with its exons and/or CDS segments
429     class GffObj:public GSeg {
430     //utility segment-merging function for addExon()
431     void expandExon(int xovl, uint segstart, uint segend,
432     char exontype, double sc, char fr, int qs, int qe);
433     protected:
434     //coordinate transformation data:
435     uint xstart; //absolute genomic coordinates of reference region
436     uint xend;
437     char xstatus; //coordinate transform status:
438     //0 : (start,end) coordinates are absolute
439     //'+' : (start,end) coords are relative to xstart..xend region
440     //'-' : (start,end) are relative to the reverse complement of xstart..xend region
441     //--
442     char* gffID; // ID name for mRNA (parent) feature
443 gpertea 16 char* gene_name; //value of gene_name attribute (GTF) if present or Name attribute of the parent gene feature (GFF3)
444     char* geneID; //value of gene_id attribute (GTF) if present or ID attribute of a parent gene feature (GFF3)
445     unsigned int flags;
446 gpertea 2 //-- friends:
447     friend class GffReader;
448     friend class GffExon;
449     public:
450 gpertea 16 static GffNames* names; // dictionary storage that holds the various attribute names etc.
451 gpertea 2 int track_id; // index of track name in names->tracks
452     int gseq_id; // index of genomic sequence name in names->gseqs
453     int ftype_id; // index of this record's feature name in names->feats, or the special gff_fid_mRNA value
454 gpertea 16 int exon_ftype_id; //index of child subfeature name in names->feats (that subfeature stored in "exons")
455 gpertea 2 //if ftype_id==gff_fid_mRNA then this value is ignored
456     GList<GffExon> exons; //for non-mRNA entries, these can be any subfeature of type subftype_id
457 gpertea 16 GPVec<GffObj> children;
458     GffObj* parent;
459 gpertea 2 int udata; //user data, flags etc.
460     void* uptr; //user pointer (to a parent object, cluster, locus etc.)
461     GffObj* ulink; //link to another GffObj (user controlled field)
462     // mRNA specific fields:
463     bool isCDS; //just a CDS, no UTRs
464     bool partial; //partial CDS
465     uint CDstart; //CDS start coord
466     uint CDend; //CDS end coord
467     char CDphase; //initial phase for CDS start
468 gpertea 16 bool hasErrors() { return ((flags & gfo_flag_HAS_ERRORS)!=0); }
469     void hasErrors(bool v) {
470     if (v) flags |= gfo_flag_HAS_ERRORS;
471     else flags &= ~gfo_flag_HAS_ERRORS;
472     }
473 gpertea 153 bool hasGffID() { return ((flags & gfo_flag_HAS_GFF_ID)!=0); }
474     void hasGffID(bool v) {
475     if (v) flags |= gfo_flag_HAS_GFF_ID;
476     else flags &= ~gfo_flag_HAS_GFF_ID;
477 gpertea 16 }
478     bool createdByExon() { return ((flags & gfo_flag_BY_EXON)!=0); }
479     void createdByExon(bool v) {
480     if (v) flags |= gfo_flag_BY_EXON;
481     else flags &= ~gfo_flag_BY_EXON;
482     }
483     bool isGene() { return ((flags & gfo_flag_IS_GENE)!=0); }
484     void isGene(bool v) {
485     if (v) flags |= gfo_flag_IS_GENE;
486     else flags &= ~gfo_flag_IS_GENE;
487     }
488     bool isDiscarded() { return ((flags & gfo_flag_DISCARDED)!=0); }
489     void isDiscarded(bool v) {
490     if (v) flags |= gfo_flag_DISCARDED;
491     else flags &= ~gfo_flag_DISCARDED;
492     }
493    
494     bool isUsed() { return ((flags & gfo_flag_LST_KEEP)!=0); }
495     void isUsed(bool v) {
496     if (v) flags |= gfo_flag_LST_KEEP;
497     else flags &= ~gfo_flag_LST_KEEP;
498     }
499     bool isTranscript() { return ((flags & gfo_flag_IS_TRANSCRIPT)!=0); }
500     void isTranscript(bool v) {
501     if (v) flags |= gfo_flag_IS_TRANSCRIPT;
502     else flags &= ~gfo_flag_IS_TRANSCRIPT;
503     }
504     bool promotedChildren() { return ((flags & gfo_flag_CHILDREN_PROMOTED)!=0); }
505     void promotedChildren(bool v) {
506     if (v) flags |= gfo_flag_CHILDREN_PROMOTED;
507     else flags &= ~gfo_flag_CHILDREN_PROMOTED;
508     }
509     void setLevel(byte v) {
510     if (v==0) flags &= ~gfo_flag_LEVEL_MSK;
511     else flags &= ~(((uint)v) << gfo_flagShift_LEVEL);
512     }
513     byte incLevel() {
514     uint v=((flags & gfo_flag_LEVEL_MSK) >> gfo_flagShift_LEVEL);
515     v++;
516     flags &= ~(v << gfo_flagShift_LEVEL);
517     return v;
518     }
519     byte getLevel() {
520     return ((byte)((flags & gfo_flag_LEVEL_MSK) >> gfo_flagShift_LEVEL));
521     }
522 gpertea 2
523 gpertea 16 bool isValidTranscript() {
524     //return (ftype_id==gff_fid_mRNA && exons.Count()>0);
525     return (isTranscript() && exons.Count()>0);
526     }
527    
528 gpertea 2
529     int addExon(uint segstart, uint segend, double sc=0, char fr='.',
530     int qs=0, int qe=0, bool iscds=false, char exontype=0);
531    
532 gpertea 16 int addExon(GffReader* reader, GffLine* gl, bool keepAttr=false, bool noExonAttr=true);
533    
534 gpertea 2 void removeExon(int idx);
535 gpertea 16 void removeExon(GffExon* p);
536 gpertea 2 char strand; //true if features are on the reverse complement strand
537     double gscore;
538     double uscore; //custom, user-computed score, if needed
539     int covlen; //total coverage of reference genomic sequence (sum of maxcf segment lengths)
540 gpertea 16
541 gpertea 2 //--------- optional data:
542     int qlen; //query length, start, end - if available
543     int qstart;
544     int qend;
545     int qcov; //query coverage - percent
546     GffAttrs* attrs; //other gff3 attributes found for the main mRNA feature
547     //constructor by gff line parsing:
548     GffObj(GffReader* gfrd, GffLine* gffline, bool keepAttrs=false, bool noExonAttr=true);
549     //if gfline->Parent!=NULL then this will also add the first sub-feature
550     // otherwise, only the main feature is created
551     void clearAttrs() {
552 gpertea 16 if (attrs!=NULL) {
553     bool sharedattrs=(exons.Count()>0 && exons[0]->attrs==attrs);
554     delete attrs; attrs=NULL;
555     if (sharedattrs) exons[0]->attrs=NULL;
556     }
557 gpertea 2 }
558 gpertea 16 GffObj(char* anid=NULL):GSeg(0,0), exons(true,true,false), children(1,false) {
559     //exons: sorted, free, non-unique
560 gpertea 2 gffID=NULL;
561     uptr=NULL;
562     ulink=NULL;
563 gpertea 16 flags=0;
564 gpertea 2 udata=0;
565 gpertea 16 parent=NULL;
566 gpertea 2 ftype_id=-1;
567 gpertea 16 exon_ftype_id=-1;
568 gpertea 2 if (anid!=NULL) gffID=Gstrdup(anid);
569     gffnames_ref(names);
570     qlen=0;
571     qstart=0;
572     qend=0;
573     qcov=0;
574     partial=true;
575     isCDS=false;
576     CDstart=0; // hasCDS <=> CDstart>0
577     CDend=0;
578     CDphase=0;
579     gseq_id=-1;
580     track_id=-1;
581     xstart=0;
582     xend=0;
583     xstatus=0;
584 gpertea 16 strand='.';
585 gpertea 2 gscore=0;
586     uscore=0;
587     attrs=NULL;
588     covlen=0;
589 gpertea 16 gene_name=NULL;
590     geneID=NULL;
591 gpertea 2 }
592     ~GffObj() {
593     GFREE(gffID);
594 gpertea 16 GFREE(gene_name);
595     GFREE(geneID);
596     clearAttrs();
597 gpertea 2 gffnames_unref(names);
598     }
599     //--------------
600 gpertea 16 GffObj* finalize(GffReader* gfr, bool mergeCloseExons=false,
601     bool keepAttrs=false, bool noExonAttr=true);
602     //complete parsing: must be called in order to merge adjacent/close proximity subfeatures
603     void parseAttrs(GffAttrs*& atrlist, char* info, bool isExon=false);
604 gpertea 2 const char* getSubfName() { //returns the generic feature type of the entries in exons array
605 gpertea 16 int sid=exon_ftype_id;
606 gpertea 2 if (sid==gff_fid_exon && isCDS) sid=gff_fid_CDS;
607     return names->feats.getName(sid);
608     }
609 gpertea 55 void addCDS(uint cd_start, uint cd_end, char phase=0);
610    
611 gpertea 16 bool monoFeature() {
612     return (exons.Count()==0 ||
613 gpertea 63 (exons.Count()==1 && //exon_ftype_id==ftype_id &&
614 gpertea 50 exons[0]->end==this->end && exons[0]->start==this->start));
615 gpertea 16 }
616    
617     bool hasCDS() { return (CDstart>0); }
618    
619 gpertea 2 const char* getFeatureName() {
620     return names->feats.getName(ftype_id);
621     }
622 gpertea 16 void setFeatureName(const char* feature);
623    
624     void addAttr(const char* attrname, const char* attrvalue);
625     int removeAttr(const char* attrname, const char* attrval=NULL);
626     int removeAttr(int aid, const char* attrval=NULL);
627     int removeExonAttr(GffExon& exon, const char* attrname, const char* attrval=NULL);
628     int removeExonAttr(GffExon& exon, int aid, const char* attrval=NULL);
629 gpertea 2 const char* getAttrName(int i) {
630     if (attrs==NULL) return NULL;
631     return names->attrs.getName(attrs->Get(i)->attr_id);
632     }
633 gpertea 16 char* getAttr(const char* attrname, bool checkFirstExon=false) {
634     if (names==NULL || attrname==NULL) return NULL;
635     char* r=NULL;
636     if (attrs==NULL) {
637     if (!checkFirstExon) return NULL;
638     }
639     else r=attrs->getAttr(names, attrname);
640     if (r!=NULL) return r;
641     if (checkFirstExon && exons.Count()>0) {
642     r=exons[0]->getAttr(names, attrname);
643     }
644     return r;
645 gpertea 2 }
646    
647 gpertea 16 char* getExonAttr(GffExon* exon, const char* attrname) {
648     if (exon==NULL || attrname==NULL) return NULL;
649     return exon->getAttr(names, attrname);
650     }
651    
652     char* getExonAttr(int exonidx, const char* attrname) {
653     if (exonidx<0 || exonidx>=exons.Count() || attrname==NULL) return NULL;
654     return exons[exonidx]->getAttr(names, attrname);
655     }
656    
657 gpertea 2 char* getAttrValue(int i) {
658     if (attrs==NULL) return NULL;
659     return attrs->Get(i)->attr_val;
660     }
661     const char* getGSeqName() {
662     return names->gseqs.getName(gseq_id);
663     }
664 gpertea 16
665     const char* getRefName() {
666     return names->gseqs.getName(gseq_id);
667     }
668     void setRefName(const char* newname);
669    
670 gpertea 2 const char* getTrackName() {
671     return names->tracks.getName(track_id);
672     }
673 gpertea 16 bool exonOverlap(uint s, uint e) {//check if ANY exon overlaps given segment
674 gpertea 2 //ignores strand!
675 gpertea 144 if (s>e) Gswap(s,e);
676 gpertea 2 for (int i=0;i<exons.Count();i++) {
677     if (exons[i]->overlap(s,e)) return true;
678     }
679     return false;
680     }
681     bool exonOverlap(GffObj& m) {//check if ANY exon overlaps given segment
682     //if (gseq_id!=m.gseq_id) return false;
683     // ignores strand and gseq_id, must check in advance
684     for (int i=0;i<exons.Count();i++) {
685     for (int j=0;j<m.exons.Count();j++) {
686     if (exons[i]->start>m.exons[j]->end) continue;
687     if (m.exons[j]->start>exons[i]->end) break;
688     //-- overlap if we are here:
689     return true;
690     }
691     }
692     return false;
693     }
694 gpertea 16
695 gpertea 2 int exonOverlapIdx(uint s, uint e, int* ovlen=NULL) {
696 gpertea 16 //return the exons' index for the overlapping OR ADJACENT exon
697 gpertea 2 //ovlen, if given, will return the overlap length
698 gpertea 144 if (s>e) Gswap(s,e);
699 gpertea 16 s--;e++; //to also catch adjacent exons
700 gpertea 2 for (int i=0;i<exons.Count();i++) {
701     if (exons[i]->start>e) break;
702     if (s>exons[i]->end) continue;
703     //-- overlap if we are here:
704     if (ovlen!=NULL) {
705 gpertea 16 s++;e--;
706 gpertea 2 int ovlend= (exons[i]->end>e) ? e : exons[i]->end;
707     *ovlen= ovlend - ((s>exons[i]->start)? s : exons[i]->start)+1;
708     }
709     return i;
710     } //for each exon
711     *ovlen=0;
712     return -1;
713     }
714 gpertea 16
715 gpertea 2 int exonOverlapLen(GffObj& m) {
716     if (start>m.end || m.start>end) return 0;
717     int i=0;
718     int j=0;
719     int ovlen=0;
720     while (i<exons.Count() && j<m.exons.Count()) {
721     uint istart=exons[i]->start;
722     uint iend=exons[i]->end;
723     uint jstart=m.exons[j]->start;
724     uint jend=m.exons[j]->end;
725     if (istart>jend) { j++; continue; }
726     if (jstart>iend) { i++; continue; }
727     //exon overlap
728     uint ovstart=GMAX(istart,jstart);
729     if (iend<jend) {
730     ovlen+=iend-ovstart+1;
731     i++;
732     }
733     else {
734     ovlen+=jend-ovstart+1;
735     j++;
736     }
737     }//while comparing exons
738     return ovlen;
739     }
740    
741     bool exonOverlap(GffObj* m) {
742     return exonOverlap(*m);
743     }
744     //---------- coordinate transformation
745     void xcoord(uint grstart, uint grend, char xstrand='+') {
746     //relative coordinate transform, and reverse-complement transform if xstrand is '-'
747     //does nothing if xstatus is the same already
748     if (xstatus) {
749     if (xstatus==xstrand && grstart==xstart && grend==xend) return;
750     unxcoord();//restore original coordinates
751     }
752     xstatus=xstrand;
753     xstart=grstart;
754     xend=grend;
755     if (CDstart>0) xcoordseg(CDstart, CDend);
756     for (int i=0;i<exons.Count();i++) {
757     xcoordseg(exons[i]->start, exons[i]->end);
758     }
759     if (xstatus=='-') {
760     exons.Reverse();
761     int flen=end-start;
762     start=xend-end+1;
763     end=start+flen;
764     }
765     else {
766     start=start-xstart+1;
767     end=end-xstart+1;
768     }
769     }
770    
771     //transform an arbitrary segment based on current xstatus/xstart-xend
772     void xcoordseg(uint& segstart, uint &segend) {
773     if (xstatus==0) return;
774     if (xstatus=='-') {
775     int flen=segend-segstart;
776     segstart=xend-segend+1;
777     segend=segstart+flen;
778     return;
779     }
780     else {
781     segstart=segstart-xstart+1;
782     segend=segend-xstart+1;
783     }
784     }
785    
786     void unxcoord() { //revert back to absolute genomic/gff coordinates if xstatus==true
787     if (xstatus==0) return; //nothing to do, no transformation appplied
788     if (CDstart>0) unxcoordseg(CDstart, CDend);
789     //restore all GffExon intervals too
790     for (int i=0;i<exons.Count();i++) {
791     unxcoordseg(exons[i]->start, exons[i]->end);
792     }
793     if (xstatus=='-') {
794     exons.Reverse();
795     int flen=end-start;
796     start=xend-end+1;
797     end=start+flen;
798     }
799     else {
800     start=start+xstart-1;
801     end=end+xstart-1;
802     }
803     xstatus=0;
804     }
805     void unxcoordseg(uint& astart, uint &aend) {
806     //restore an arbitrary interval -- does NOT change the transform state!
807     if (xstatus==0) return;
808     if (xstatus=='-') {
809     int flen=aend-astart;
810     astart=xend-aend+1;
811     aend=astart+flen;
812     }
813     else {
814     astart=astart+xstart-1;
815     aend=aend+xstart-1;
816     }
817     }
818     //---------------------
819     bool operator==(GffObj& d){
820     return (gseq_id==d.gseq_id && start==d.start && end==d.end && strcmp(gffID, d.gffID)==0);
821     }
822     bool operator>(GffObj& d){
823     if (gseq_id!=d.gseq_id) return (gseq_id>d.gseq_id);
824     if (start==d.start) {
825 gpertea 16 if (getLevel()==d.getLevel()) {
826     if (end==d.end) return (strcmp(gffID, d.gffID)>0);
827     else return (end>d.end);
828     } else return (getLevel()>d.getLevel());
829 gpertea 2 } else return (start>d.start);
830     }
831     bool operator<(GffObj& d){
832     if (gseq_id!=d.gseq_id) return (gseq_id<d.gseq_id);
833     if (start==d.start) {
834 gpertea 16 if (getLevel()==d.getLevel()) {
835     if (end==d.end) return strcmp(gffID, d.gffID)<0;
836 gpertea 2 else return end<d.end;
837 gpertea 16 } else return (getLevel()<d.getLevel());
838 gpertea 2 } else return (start<d.start);
839     }
840     char* getID() { return gffID; }
841 gpertea 16 char* getGeneID() { return geneID; }
842     char* getGeneName() { return gene_name; }
843     void setGeneName(const char* gname) {
844     GFREE(gene_name);
845     if (gname) gene_name=Gstrdup(gname);
846     }
847     void setGeneID(const char* gene_id) {
848     GFREE(geneID);
849     if (gene_id) geneID=Gstrdup(gene_id);
850     }
851 gpertea 2 int addSeg(GffLine* gfline);
852     int addSeg(int fnid, GffLine* gfline);
853     void getCDSegs(GArray<GffCDSeg>& cds);
854 gpertea 16
855     void updateExonPhase(); //for CDS-only features, updates GExon::phase
856    
857     void printGxfLine(FILE* fout, const char* tlabel, const char* gseqname,
858     bool iscds, uint segstart, uint segend, int exidx, char phase, bool gff3);
859     void printGxf(FILE* fout, GffPrintMode gffp=pgffExon,
860     const char* tlabel=NULL, const char* gfparent=NULL);
861     void printGtf(FILE* fout, const char* tlabel=NULL) {
862 gpertea 2 printGxf(fout, pgtfAny, tlabel);
863     }
864 gpertea 16 void printGff(FILE* fout, const char* tlabel=NULL,
865     const char* gfparent=NULL) {
866     printGxf(fout, pgffAny, tlabel, gfparent);
867 gpertea 2 }
868 gpertea 16 void printTranscriptGff(FILE* fout, char* tlabel=NULL,
869     bool showCDS=false, const char* gfparent=NULL) {
870     if (isValidTranscript())
871     printGxf(fout, showCDS ? pgffBoth : pgffExon, tlabel, gfparent);
872 gpertea 2 }
873     void printSummary(FILE* fout=NULL);
874     void getCDS_ends(uint& cds_start, uint& cds_end);
875     void mRNA_CDS_coords(uint& cds_start, uint& cds_end);
876     char* getSpliced(GFaSeqGet* faseq, bool CDSonly=false, int* rlen=NULL,
877     uint* cds_start=NULL, uint* cds_end=NULL, GList<GSeg>* seglst=NULL);
878 gpertea 16 char* getUnspliced(GFaSeqGet* faseq, int* rlen, GList<GSeg>* seglst);
879 gpertea 2 char* getSplicedTr(GFaSeqGet* faseq, bool CDSonly=true, int* rlen=NULL);
880     //bool validCDS(GFaSeqGet* faseq); //has In-Frame Stop Codon ?
881     bool empty() { return (start==0); }
882     };
883    
884     typedef bool GffRecFunc(GffObj* gobj, void* usrptr1, void* usrptr2);
885     //user callback after parsing a mapping object:
886     // Returns: "done with it" status:
887     // TRUE if gobj is no longer needed so it's FREEd upon return
888     // FALSE if the user needs the gobj pointer and is responsible for
889     // collecting and freeing all GffObj objects
890    
891    
892     //GSeqStat: collect basic stats about a common underlying genomic sequence
893     // for multiple GffObj
894     class GSeqStat {
895     public:
896     int gseqid; //gseq id in the global static pool of gseqs
897     char* gseqname; //just a pointer to the name of gseq
898     //int fcount;//number of features on this gseq
899     uint mincoord;
900     uint maxcoord;
901 gpertea 16 uint maxfeat_len; //maximum feature length on this genomic sequence
902     GffObj* maxfeat;
903     GSeqStat(int id=-1, char* name=NULL) {
904 gpertea 2 gseqid=id;
905     gseqname=name;
906     mincoord=MAXUINT;
907     maxcoord=0;
908 gpertea 16 maxfeat_len=0;
909     maxfeat=NULL;
910 gpertea 2 }
911     bool operator>(GSeqStat& g) {
912     return (gseqid>g.gseqid);
913     }
914     bool operator<(GSeqStat& g) {
915     return (gseqid<g.gseqid);
916     }
917     bool operator==(GSeqStat& g) {
918     return (gseqid==g.gseqid);
919     }
920     };
921    
922    
923     int gfo_cmpByLoc(const pointer p1, const pointer p2);
924    
925     class GfList: public GList<GffObj> {
926 gpertea 16 //just adding the option to sort by genomic sequence and coordinate
927     bool mustSort;
928 gpertea 2 public:
929     GfList(bool sortbyloc=false):GList<GffObj>(false,false,false) {
930 gpertea 16 //GffObjs in this list are NOT deleted when the list is cleared
931     //-- for deallocation of these objects, call freeAll() or freeUnused() as needed
932     mustSort=sortbyloc;
933     }
934     void sortedByLoc(bool v=true) {
935     bool prev=mustSort;
936     mustSort=v;
937     if (fCount>0 && mustSort && !prev) {
938     this->setSorted((GCompareProc*)gfo_cmpByLoc);
939     }
940     }
941     void finalize(GffReader* gfr, bool mergeCloseExons,
942     bool keepAttrs=false, bool noExonAttr=true) { //if set, enforce sort by locus
943     if (mustSort) { //force (re-)sorting
944     this->setSorted(false);
945     this->setSorted((GCompareProc*)gfo_cmpByLoc);
946     }
947     int delcount=0;
948     for (int i=0;i<Count();i++) {
949     //finish the parsing for each GffObj
950     fList[i]->finalize(gfr, mergeCloseExons, keepAttrs, noExonAttr);
951     }
952     if (delcount>0) this->Pack();
953     }
954     void freeAll() {
955     for (int i=0;i<fCount;i++) {
956     delete fList[i];
957     fList[i]=NULL;
958     }
959     Clear();
960     }
961     void freeUnused() {
962     for (int i=0;i<fCount;i++) {
963     if (fList[i]->isUsed()) continue;
964     //inform the children
965     for (int c=0;c<fList[i]->children.Count();c++) {
966     fList[i]->children[c]->parent=NULL;
967     }
968     delete fList[i];
969     fList[i]=NULL;
970     }
971     Clear();
972     }
973    
974 gpertea 2 };
975    
976 gpertea 151 struct GfoHolder {
977 gpertea 16 int idx; //position in GffReader::gflst array
978 gpertea 151 GffObj* gffobj;
979 gpertea 16 GfoHolder(GffObj* gfo=NULL, int i=0) {
980     idx=i;
981     gffobj=gfo;
982     }
983     };
984    
985     class CNonExon { //utility class used in subfeature promotion
986     public:
987     int idx;
988     GffObj* parent;
989     GffExon* exon;
990     GffLine* gffline;
991     CNonExon(int i, GffObj* p, GffExon* e, GffLine* gl) {
992     parent=p;
993     exon=e;
994     idx=i;
995     gffline=new GffLine(gl);
996     }
997     ~CNonExon() {
998     delete gffline;
999     }
1000     };
1001    
1002    
1003 gpertea 2 class GffReader {
1004     friend class GffObj;
1005     friend class GffLine;
1006     char* linebuf;
1007     off_t fpos;
1008     int buflen;
1009     protected:
1010 gpertea 16 bool gff_warns; //warn about duplicate IDs, etc. even when they are on different chromosomes
1011 gpertea 2 FILE* fh;
1012     char* fname; //optional fasta file with the underlying genomic sequence to be attached to this reader
1013     GffNames* names; //just a pointer to the global static Gff names repository in GffObj
1014     GffLine* gffline;
1015 gpertea 16 bool transcriptsOnly; //keep only transcripts w/ their exon/CDS features
1016     GHash<int> discarded_ids; //for transcriptsOnly mode, keep track
1017     // of discarded parent IDs
1018 gpertea 151 GHash< GVec<GfoHolder> > phash; //transcript_id+contig (Parent~Contig) => [gflst index, GffObj]
1019 gpertea 153 //GHash<int> tids; //just for transcript_id uniqueness
1020 gpertea 2 char* gfoBuildId(const char* id, const char* ctg);
1021 gpertea 151 //void gfoRemove(const char* id, const char* ctg);
1022     GfoHolder* gfoAdd(GffObj* gfo, int idx);
1023     GfoHolder* gfoAdd(GVec<GfoHolder>& glst, GffObj* gfo, int idx);
1024 gpertea 195 // const char* id, const char* ctg, char strand, GVec<GfoHolder>** glst, uint start, uint end
1025     GfoHolder* gfoFind(const char* id, const char* ctg=NULL, GVec<GfoHolder>** glst=NULL,
1026     char strand=0, uint start=0, uint end=0);
1027 gpertea 16 CNonExon* subfPoolCheck(GffLine* gffline, GHash<CNonExon>& pex, char*& subp_name);
1028     void subfPoolAdd(GHash<CNonExon>& pex, GfoHolder* newgfo);
1029     GfoHolder* promoteFeature(CNonExon* subp, char*& subp_name, GHash<CNonExon>& pex,
1030     bool keepAttr, bool noExonAttr);
1031 gpertea 2 public:
1032 gpertea 16 GfList gflst; //accumulate GffObjs being read
1033     GfoHolder* newGffRec(GffLine* gffline, bool keepAttr, bool noExonAttr,
1034 gpertea 151 GffObj* parent=NULL, GffExon* pexon=NULL, GVec<GfoHolder>* glst=NULL);
1035 gpertea 16 GfoHolder* replaceGffRec(GffLine* gffline, bool keepAttr, bool noExonAttr, int replaceidx);
1036     GfoHolder* updateGffRec(GfoHolder* prevgfo, GffLine* gffline,
1037     bool keepAttr);
1038     GfoHolder* updateParent(GfoHolder* newgfh, GffObj* parent);
1039     bool addExonFeature(GfoHolder* prevgfo, GffLine* gffline, GHash<CNonExon>& pex, bool noExonAttr);
1040 gpertea 2 GList<GSeqStat> gseqstats; //list of all genomic sequences seen by this reader, accumulates stats
1041 gpertea 16 GffReader(FILE* f=NULL, bool t_only=false, bool sortbyloc=false):discarded_ids(true),
1042 gpertea 153 phash(true), gflst(sortbyloc), gseqstats(true,true,true) {
1043 gpertea 16 gff_warns=gff_show_warnings;
1044 gpertea 2 names=NULL;
1045     gffline=NULL;
1046 gpertea 16 transcriptsOnly=t_only;
1047 gpertea 2 fpos=0;
1048     fname=NULL;
1049     fh=f;
1050     GMALLOC(linebuf, GFF_LINELEN);
1051     buflen=GFF_LINELEN-1;
1052     }
1053 gpertea 16 void init(FILE *f, bool t_only=false, bool sortbyloc=false) {
1054     fname=NULL;
1055     fh=f;
1056     if (fh!=NULL) rewind(fh);
1057     fpos=0;
1058     transcriptsOnly=t_only;
1059     gflst.sortedByLoc(sortbyloc);
1060     }
1061     GffReader(char* fn, bool t_only=false, bool sort=false):discarded_ids(true), phash(true),
1062 gpertea 153 gflst(sort),gseqstats(true,true,true) {
1063 gpertea 16 gff_warns=gff_show_warnings;
1064 gpertea 2 names=NULL;
1065     fname=Gstrdup(fn);
1066 gpertea 16 transcriptsOnly=t_only;
1067 gpertea 2 fh=fopen(fname, "rb");
1068     fpos=0;
1069     gffline=NULL;
1070     GMALLOC(linebuf, GFF_LINELEN);
1071     buflen=GFF_LINELEN-1;
1072     }
1073    
1074 gpertea 16 ~GffReader() {
1075     delete gffline;
1076 gpertea 2 gffline=NULL;
1077     fpos=0;
1078 gpertea 16 gflst.freeUnused();
1079 gpertea 2 gflst.Clear();
1080 gpertea 16 discarded_ids.Clear();
1081 gpertea 2 phash.Clear();
1082     gseqstats.Clear();
1083     GFREE(fname);
1084     GFREE(linebuf);
1085     }
1086    
1087 gpertea 16 void showWarnings(bool v=true) {
1088     gff_warns=v;
1089     gff_show_warnings=v;
1090     }
1091    
1092 gpertea 2 GffLine* nextGffLine();
1093    
1094 gpertea 16 // load all subfeatures, re-group them:
1095     void readAll(bool keepAttr=false, bool mergeCloseExons=false, bool noExonAttr=true);
1096 gpertea 2
1097     }; // end of GffReader
1098    
1099     #endif