1 |
+ |
#include "gff_utils.h" |
2 |
|
#include "GArgs.h" |
2 |
– |
#include "GStr.h" |
3 |
– |
#include "GHash.hh" |
4 |
– |
#include "GList.hh" |
5 |
– |
#include "GFaSeqGet.h" |
6 |
– |
#include "gff.h" |
3 |
|
#include <ctype.h> |
4 |
< |
// do not care about cdb compression |
5 |
< |
#ifdef ENABLE_COMPRESSION |
6 |
< |
#undef ENABLE_COMPRESSION |
7 |
< |
#endif |
8 |
< |
#include "GCdbYank.h" |
4 |
> |
// don't care about cdb compression |
5 |
> |
//#ifdef ENABLE_COMPRESSION |
6 |
> |
//#undef ENABLE_COMPRESSION |
7 |
> |
//#endif |
8 |
> |
//#include "GCdbYank.h" |
9 |
|
|
10 |
|
#define USAGE "Usage:\n\ |
11 |
< |
gffread <input_gff..> [-g <genomic_seq_fasta> | <dir>][-s <seq_info.fsize>]\n\ |
12 |
< |
[-o <outfile.gff3>] [-t <tname>] [-r [<chr>:]<start>..<end>] [-i <maxintron>]\n\ |
13 |
< |
[-CTVNMAFGRUVBHSZWTO] [-w <spl_exons.fa>] [-x <spl_cds.fa>] [-y <tran_cds.fa>]\n\ |
14 |
< |
\n\ |
15 |
< |
Filters and/or converts GFF/GTF records.\n\ |
11 |
> |
gffread <input_gff> [-g <genomic_seqs_fasta> | <dir>][-s <seq_info.fsize>] \n\ |
12 |
> |
[-o <outfile.gff>] [-t <tname>] [-r [[<strand>]<chr>:]<start>..<end>] \n\ |
13 |
> |
[-CTVNJMKQAFGRUVBHZWTOLE] [-w <spl_exons.fa>] [-x <spl_cds.fa>] [-y <tr_cds.fa>]\n\ |
14 |
> |
[-i <maxintron>] \n\ |
15 |
> |
Filters and/or converts GFF3/GTF2 records.\n\ |
16 |
> |
<input_gff> is a GFF file, use '-' if the GFF records will be given at stdin\n\ |
17 |
|
\n\ |
18 |
|
Options:\n\ |
19 |
< |
-g full path to one fasta file with the genomic sequence\n\ |
20 |
< |
for all input mappings, OR a path to the directory where\n\ |
21 |
< |
genomic sequences can be found as single-fasta files\n\ |
22 |
< |
-s for mRNA/EST/protein mappings a tab-delimited file provides this info\n\ |
19 |
> |
-g full path to a multi-fasta file with the genomic sequences\n\ |
20 |
> |
for all input mappings, OR a directory with single-fasta files\n\ |
21 |
> |
(one per genomic sequence, with file names matching sequence names)\n\ |
22 |
> |
-s <seq_info.fsize> is a tab-delimited file providing this info\n\ |
23 |
|
for each of the mapped sequences:\n\ |
24 |
|
<seq-name> <seq-length> <seq-description>\n\ |
25 |
+ |
(useful for -A option with mRNA/EST/protein mappings)\n\ |
26 |
|
-i discard transcripts having an intron larger than <maxintron>\n\ |
27 |
|
-r only show transcripts crossing coordinate range <start>..<end>\n\ |
28 |
< |
(on chromosome/contig <chr> if provided)\n\ |
28 |
> |
(on chromosome/contig <chr>, strand <strand> if provided)\n\ |
29 |
|
-R for -r option, discard all transcripts that are not fully \n\ |
30 |
|
contained within given range\n\ |
31 |
|
-U discard single-exon transcripts\n\ |
32 |
< |
-C discard mRNAs that have no CDS feature\n\ |
33 |
< |
-F keep all attributes from last column of GFF/GTF\n\ |
32 |
> |
-C coding only: discard mRNAs that have no CDS feature\n\ |
33 |
> |
-F full GFF attribute preservation (all attributes are shown)\n\ |
34 |
|
-G only parse additional exon attributes from the first exon\n\ |
35 |
|
and move them to the mRNA level (useful for GTF input)\n\ |
36 |
< |
-A use the description field from <seq_info.fsize> and add it as\n\ |
37 |
< |
a descr=.. attribute to the mRNA and/or Gene record\n\ |
36 |
> |
-A use the description field from <seq_info.fsize> and add it\n\ |
37 |
> |
as the value for a 'descr' attribute to the GFF record\n\ |
38 |
|
\n\ |
39 |
< |
-O process other (non-mRNA) GFF/GTF records (default is discard non-mRNA\n\ |
40 |
< |
data). \n\ |
43 |
< |
Note: non-mRNA records must have only one subfeature per parent feature\n\ |
39 |
> |
-O process also non-transcript GFF records (by default non-transcript\n\ |
40 |
> |
records are ignored)\n\ |
41 |
|
-V discard any mRNAs with CDS having in-frame stop codons\n\ |
42 |
|
-H for -V option, check and adjust the starting CDS phase\n\ |
43 |
|
if the original phase leads to a translation with an \n\ |
44 |
|
in-frame stop codon\n\ |
45 |
|
-B for -V option, single-exon transcripts are also checked on the\n\ |
46 |
|
opposite strand\n\ |
47 |
< |
-N only show multi-exon mRNAs if all their introns have the \n\ |
48 |
< |
typical splice site consensus ( GT-AG, GC-AG or AT-AC )\n\ |
49 |
< |
-M discard any mRNAs that either lack initial START codon\n\ |
47 |
> |
-N discard multi-exon mRNAs that have any intron with a non-canonical\n\ |
48 |
> |
splice site consensus (i.e. not GT-AG, GC-AG or AT-AC)\n\ |
49 |
> |
-J discard any mRNAs that either lack initial START codon\n\ |
50 |
|
or the terminal STOP codon, or have an in-frame stop codon\n\ |
51 |
|
(only print mRNAs with a fulll, valid CDS)\n\ |
52 |
|
\n\ |
53 |
< |
-S sort output GFF records by genomic sequence and start coordinate\n\ |
54 |
< |
note that this option is automatically enabled by -g option\n\ |
55 |
< |
-Z merge close exons into a single exon (intron size<4)\n\ |
56 |
< |
-t use <trackname> in the second column of each GFF/GTF output line\n\ |
57 |
< |
-w write a fasta file with spliced exons for each mapping\n\ |
58 |
< |
-x write a fasta file with spliced CDS for each mapping\n\ |
59 |
< |
-W for -w and -x options, also write for each fasta record the exon \n\ |
53 |
> |
-M/--merge : cluster the input transcripts into loci, collapsing matching\n\ |
54 |
> |
transcripts (those with the same exact introns and fully contained)\n\ |
55 |
> |
--cluster-only: same as --merge but without collapsing matching transcripts\n\ |
56 |
> |
-K for -M option: also collapse shorter, fully contained transcripts\n\ |
57 |
> |
with fewer introns than the container\n\ |
58 |
> |
-Q for -M option, remove the containment restriction:\n\ |
59 |
> |
(multi-exon transcripts will be collapsed if just their introns match,\n\ |
60 |
> |
while single-exon transcripts can partially overlap (80%))\n\ |
61 |
> |
\n\ |
62 |
> |
-E expose (warn about) duplicate transcript IDs and other potential \n\ |
63 |
> |
problems with the given GFF/GTF records\n\ |
64 |
> |
-Z merge close exons into a single exon (for intron size<4)\n\ |
65 |
> |
-w write a fasta file with spliced exons for each GFF transcript\n\ |
66 |
> |
-x write a fasta file with spliced CDS for each GFF transcript\n\ |
67 |
> |
-W for -w and -x options, also write for each fasta record the exon\n\ |
68 |
|
coordinates projected onto the spliced sequence\n\ |
69 |
< |
-y write a protein fasta file with the CDS translation for each mapping\n\ |
70 |
< |
-o the output gff3 file with the 'filtered' entries\n\ |
71 |
< |
-T output GTF format instead of GFF3\n\ |
72 |
< |
\n\ |
73 |
< |
" |
69 |
> |
-y write a protein fasta file with the translation of CDS for each record\n\ |
70 |
> |
-L Ensembl GTF to GFF3 conversion (implies -F; should be used with -m)\n\ |
71 |
> |
-m <chr_replace> is a reference (genomic) sequence replacement table with\n\ |
72 |
> |
this format:\n\ |
73 |
> |
<original_ref_ID> <new_ref_ID>\n\ |
74 |
> |
GFF records on reference sequences that are not found among the\n\ |
75 |
> |
<original_ref_ID> entries in this file will be filtered out\n\ |
76 |
> |
-o the \"filtered\" GFF records will be written to <outfile.gff>\n\ |
77 |
> |
(use -o- for printing to stdout)\n\ |
78 |
> |
-t use <trackname> in the second column of each GFF output line\n\ |
79 |
> |
-T -o option will output GTF format instead of GFF3\n\ |
80 |
> |
" |
81 |
> |
|
82 |
> |
|
83 |
> |
class SeqInfo { //populated from the -s option of gffread |
84 |
> |
public: |
85 |
> |
int len; |
86 |
> |
char* descr; |
87 |
> |
SeqInfo( int l, char* s) { |
88 |
> |
len=l; |
89 |
> |
if (s==NULL) { |
90 |
> |
descr=NULL; |
91 |
> |
} else { |
92 |
> |
descr=Gstrdup(s); |
93 |
> |
} |
94 |
> |
} |
95 |
> |
~SeqInfo() { |
96 |
> |
GFREE(descr); |
97 |
> |
} |
98 |
> |
}; |
99 |
> |
|
100 |
> |
class RefTran { |
101 |
> |
public: |
102 |
> |
char* new_name; |
103 |
> |
RefTran(char *ns) { |
104 |
> |
new_name=NULL; |
105 |
> |
if (ns!=NULL) |
106 |
> |
new_name=Gstrdup(ns); |
107 |
> |
} |
108 |
> |
~RefTran() { |
109 |
> |
GFREE(new_name); |
110 |
> |
} |
111 |
> |
}; |
112 |
|
|
113 |
|
FILE* ffasta=NULL; |
114 |
|
FILE* f_in=NULL; |
127 |
|
|
128 |
|
bool fullCDSonly=false; // starts with START, ends with STOP codon |
129 |
|
bool fullattr=false; |
130 |
< |
bool sortByLoc=false; // if the GFF output should be sorted by location |
131 |
< |
GStr gseqpath; |
132 |
< |
bool multiGSeq=false; //if a directory or a .cidx file was given to -g option |
130 |
> |
//bool sortByLoc=false; // if the GFF output should be sorted by location |
131 |
> |
bool ensembl_convert=false; //-L, assisst in converting Ensembl GTF to GFF3 |
132 |
> |
|
133 |
> |
|
134 |
> |
//GStr gseqpath; |
135 |
> |
//GStr gcdbfa; |
136 |
> |
//bool multiGSeq=false; //if a directory or a .cidx file was given to -g option |
137 |
> |
//GFaSeqGet* faseq=NULL; |
138 |
> |
//GCdbYank* gcdb=NULL; |
139 |
> |
//int gseq_id=-1; //current genome sequence ID -- the current GffObj::gseq_id |
140 |
|
bool fmtGTF=false; |
91 |
– |
int gseq_id=-1; //current genome sequence ID -- the current GffObj::gseq_id |
92 |
– |
GFaSeqGet* faseq=NULL; |
93 |
– |
GCdbYank* gcdb=NULL; |
94 |
– |
GStr gcdbfa; |
141 |
|
bool addDescr=false; |
142 |
|
//bool protmap=false; |
143 |
|
bool multiExon=false; |
147 |
|
bool mergeCloseExons=false; |
148 |
|
//range filter: |
149 |
|
char* rfltGSeq=NULL; |
150 |
+ |
char rfltStrand=0; |
151 |
|
uint rfltStart=0; |
152 |
|
uint rfltEnd=MAX_UINT; |
153 |
|
bool rfltWithin=false; //check for full containment within given range |
154 |
|
bool noExonAttr=false; |
108 |
– |
class SeqInfo { |
109 |
– |
public: |
110 |
– |
int len; |
111 |
– |
char* descr; |
112 |
– |
SeqInfo( int l, char* s) { |
113 |
– |
len=l; |
114 |
– |
if (s==NULL) { |
115 |
– |
descr=NULL; |
116 |
– |
} else { |
117 |
– |
descr=Gstrdup(s); |
118 |
– |
} |
119 |
– |
} |
120 |
– |
~SeqInfo() { |
121 |
– |
GFREE(descr); |
122 |
– |
} |
123 |
– |
}; |
155 |
|
|
156 |
+ |
bool doCluster=false; |
157 |
+ |
bool doCollapseRedundant=false; |
158 |
|
|
159 |
< |
class GSpliceSite { |
127 |
< |
public: |
128 |
< |
char nt[3]; |
129 |
< |
GSpliceSite(const char* c, bool revc=false) { |
130 |
< |
nt[2]=0; |
131 |
< |
if (c==NULL) { |
132 |
< |
nt[0]=0; |
133 |
< |
nt[1]=0; |
134 |
< |
return; |
135 |
< |
} |
136 |
< |
if (revc) { |
137 |
< |
nt[0]=toupper(ntComplement(c[1])); |
138 |
< |
nt[1]=toupper(ntComplement(c[0])); |
139 |
< |
} |
140 |
< |
else { |
141 |
< |
nt[0]=toupper(c[0]); |
142 |
< |
nt[1]=toupper(c[1]); |
143 |
< |
} |
144 |
< |
} |
145 |
< |
|
146 |
< |
GSpliceSite(const char* intron, int intronlen, bool getAcceptor, bool revc=false) { |
147 |
< |
nt[2]=0; |
148 |
< |
if (intron==NULL || intronlen==0) |
149 |
< |
GError("Error: invalid intron or intron len for GSpliceSite()!\n"); |
150 |
< |
const char* c=intron; |
151 |
< |
if (revc) { |
152 |
< |
if (!getAcceptor) c+=intronlen-2; |
153 |
< |
nt[0]=toupper(ntComplement(c[1])); |
154 |
< |
nt[1]=toupper(ntComplement(c[0])); |
155 |
< |
} |
156 |
< |
else { //on forward strand |
157 |
< |
if (getAcceptor) c+=intronlen-2; |
158 |
< |
nt[0]=toupper(c[0]); |
159 |
< |
nt[1]=toupper(c[1]); |
160 |
< |
}//forward strand |
161 |
< |
} |
162 |
< |
|
163 |
< |
GSpliceSite(const char n1, const char n2) { |
164 |
< |
nt[2]=0; |
165 |
< |
nt[0]=toupper(n1); |
166 |
< |
nt[1]=toupper(n2); |
167 |
< |
} |
168 |
< |
bool canonicalDonor() { |
169 |
< |
return (nt[0]=='G' && (nt[1]=='C' || nt[1]=='T')); |
170 |
< |
} |
171 |
< |
bool operator==(GSpliceSite& c) { |
172 |
< |
return (c.nt[0]==nt[0] && c.nt[1]==nt[1]); |
173 |
< |
} |
174 |
< |
bool operator==(GSpliceSite* c) { |
175 |
< |
return (c->nt[0]==nt[0] && c->nt[1]==nt[1]); |
176 |
< |
} |
177 |
< |
bool operator==(const char* c) { |
178 |
< |
//return (nt[0]==toupper(c[0]) && nt[1]==toupper(c[1])); |
179 |
< |
//assumes given const nucleotides are uppercase already! |
180 |
< |
return (nt[0]==c[0] && nt[1]==c[1]); |
181 |
< |
} |
182 |
< |
bool operator!=(const char* c) { |
183 |
< |
//assumes given const nucleotides are uppercase already! |
184 |
< |
return (nt[0]!=c[0] || nt[1]!=c[1]); |
185 |
< |
} |
186 |
< |
}; |
159 |
> |
GList<GenomicSeqData> g_data(true,true,true); //list of GFF records by genomic seq |
160 |
|
|
161 |
|
//hash with sequence info |
162 |
|
GHash<SeqInfo> seqinfo; |
163 |
|
GHash<int> isoCounter; //counts the valid isoforms |
164 |
+ |
GHash<RefTran> reftbl; |
165 |
+ |
GHash<GeneInfo> gene_ids; |
166 |
+ |
//min-max gene span associated to chr|gene_id (mostly for Ensembl conversion) |
167 |
|
|
168 |
|
bool debugMode=false; |
169 |
|
bool verbose=false; |
170 |
|
|
171 |
+ |
void loadSeqInfo(FILE* f, GHash<SeqInfo> &si) { |
172 |
+ |
GLineReader fr(f); |
173 |
+ |
while (!fr.isEof()) { |
174 |
+ |
char* line=fr.getLine(); |
175 |
+ |
if (line==NULL) break; |
176 |
+ |
char* id=line; |
177 |
+ |
char* lenstr=NULL; |
178 |
+ |
char* text=NULL; |
179 |
+ |
char* p=line; |
180 |
+ |
while (*p!=0 && !isspace(*p)) p++; |
181 |
+ |
if (*p==0) continue; |
182 |
+ |
*p=0;p++; |
183 |
+ |
while (*p==' ' || *p=='\t') p++; |
184 |
+ |
if (*p==0) continue; |
185 |
+ |
lenstr=p; |
186 |
+ |
while (*p!=0 && !isspace(*p)) p++; |
187 |
+ |
if (*p!=0) { *p=0;p++; } |
188 |
+ |
while (*p==' ' || *p=='\t') p++; |
189 |
+ |
if (*p!=0) text=p; //else text remains NULL |
190 |
+ |
int len=0; |
191 |
+ |
if (!parseInt(lenstr,len)) { |
192 |
+ |
GMessage("Warning: could not parse sequence length: %s %s\n", |
193 |
+ |
id, lenstr); |
194 |
+ |
continue; |
195 |
+ |
} |
196 |
+ |
// --- here we have finished parsing the line |
197 |
+ |
si.Add(id, new SeqInfo(len,text)); |
198 |
+ |
} //while lines |
199 |
+ |
} |
200 |
+ |
|
201 |
+ |
void loadRefTable(FILE* f, GHash<RefTran>& rt) { |
202 |
+ |
GLineReader fr(f); |
203 |
+ |
char* line=NULL; |
204 |
+ |
while ((line=fr.getLine())) { |
205 |
+ |
char* orig_id=line; |
206 |
+ |
char* p=line; |
207 |
+ |
while (*p!=0 && !isspace(*p)) p++; |
208 |
+ |
if (*p==0) continue; |
209 |
+ |
*p=0;p++;//split the line here |
210 |
+ |
while (*p==' ' || *p=='\t') p++; |
211 |
+ |
if (*p==0) continue; |
212 |
+ |
rt.Add(orig_id, new RefTran(p)); |
213 |
+ |
} //while lines |
214 |
+ |
} |
215 |
+ |
|
216 |
|
char* getSeqDescr(char* seqid) { |
217 |
|
static char charbuf[128]; |
218 |
|
if (seqinfo.Count()==0) return NULL; |
252 |
|
return charbuf; |
253 |
|
} |
254 |
|
|
255 |
< |
void printFasta(FILE* f, GStr& defline, char* seq, int seqlen=-1) { |
256 |
< |
if (seq==NULL) return; |
257 |
< |
int len=(seqlen>0)?seqlen:strlen(seq); |
237 |
< |
if (len<=0) return; |
238 |
< |
if (!defline.is_empty()) |
239 |
< |
fprintf(f, ">%s\n",defline.chars()); |
240 |
< |
int ilen=0; |
241 |
< |
for (int i=0; i < len; i++, ilen++) { |
242 |
< |
if (ilen == 70) { |
243 |
< |
fputc('\n', f); |
244 |
< |
ilen = 0; |
245 |
< |
} |
246 |
< |
putc(seq[i], f); |
247 |
< |
} //for |
248 |
< |
fputc('\n', f); |
249 |
< |
} |
250 |
< |
|
251 |
< |
void loadSeqInfo(FILE* f, GHash<SeqInfo> &si) { |
252 |
< |
GLineReader fr(f); |
253 |
< |
while (!fr.isEof()) { |
254 |
< |
char* line=fr.getLine(); |
255 |
< |
if (line==NULL) break; |
256 |
< |
char* id=line; |
257 |
< |
char* lenstr=NULL; |
258 |
< |
char* text=NULL; |
259 |
< |
char* p=line; |
260 |
< |
while (*p!=0 && !isspace(*p)) p++; |
261 |
< |
if (*p==0) continue; |
262 |
< |
*p=0;p++; |
263 |
< |
while (*p==' ' || *p=='\t') p++; |
264 |
< |
if (*p==0) continue; |
265 |
< |
lenstr=p; |
266 |
< |
while (*p!=0 && !isspace(*p)) p++; |
267 |
< |
if (*p!=0) { *p=0;p++; } |
268 |
< |
while (*p==' ' || *p=='\t') p++; |
269 |
< |
if (*p!=0) text=p; //else text remains NULL |
270 |
< |
int len=0; |
271 |
< |
if (!parseInt(lenstr,len)) { |
272 |
< |
GMessage("Warning: could not parse sequence length: %s %s\n", |
273 |
< |
id, lenstr); |
274 |
< |
continue; |
275 |
< |
} |
276 |
< |
// --- here we have finished parsing the line |
277 |
< |
si.Add(id, new SeqInfo(len,text)); |
278 |
< |
} //while lines |
255 |
> |
GFaSeqGet* fastaSeqGet(GFastaDb& gfasta, GffObj& gffrec) { |
256 |
> |
if (gfasta.fastaPath==NULL) return NULL; |
257 |
> |
return gfasta.fetch(gffrec.gseq_id); |
258 |
|
} |
259 |
|
|
260 |
|
|
261 |
< |
GFaSeqGet* openFastaSeq(const char* gseqname) { |
283 |
< |
//this is only called when either we have a cidx or a directory in gseqpath |
284 |
< |
GFaSeqGet* r=NULL; |
285 |
< |
if (gcdb!=NULL) { |
286 |
< |
off_t rpos=gcdb->getRecordPos(gseqname); |
287 |
< |
if (rpos<0) // genomic sequence not found |
288 |
< |
GError("Error: cannot find genomic sequence '%s' in %s\n",gseqname, gseqpath.chars()); |
289 |
< |
return new GFaSeqGet(gcdbfa.chars(),rpos, false);// WARNING: does not validate FASTA line-len uniformity! |
290 |
< |
} |
291 |
< |
else { |
292 |
< |
GStr s(gseqpath); |
293 |
< |
s.append('/'); |
294 |
< |
s.append(gseqname); |
295 |
< |
if (debugMode) GMessage("[D] >> opening fasta file for genomic sequence %s..\n",gseqname); |
296 |
< |
if (fileExists(s.chars())==2) { |
297 |
< |
r=new GFaSeqGet(s.chars(),false); // WARNING: does not validate FASTA line-len uniformity! |
298 |
< |
r->loadall(); |
299 |
< |
if (debugMode) GMessage("[D] << loaded.\n"); |
300 |
< |
return r; |
301 |
< |
} |
302 |
< |
s.append(".fa"); |
303 |
< |
if (fileExists(s.chars())==2) { |
304 |
< |
r=new GFaSeqGet(s.chars(), false); // WARNING: does not validate FASTA line-len uniformity! |
305 |
< |
r->loadall(); |
306 |
< |
if (debugMode) GMessage("[D] << loaded.\n"); |
307 |
< |
return r; |
308 |
< |
} |
309 |
< |
s.append("sta"); |
310 |
< |
if (fileExists(s.chars())==2) { |
311 |
< |
return new GFaSeqGet(s.chars(), false); // WARNING: does not validate FASTA line-len uniformity! |
312 |
< |
r->loadall(); |
313 |
< |
if (debugMode) GMessage("[D] << loaded.\n"); |
314 |
< |
return r; |
315 |
< |
} |
316 |
< |
GError("Error: cannot load genomic sequence %s from directory %s !\n", |
317 |
< |
gseqname, gseqpath.chars()); |
318 |
< |
} |
319 |
< |
return NULL; |
320 |
< |
} |
321 |
< |
|
322 |
< |
|
323 |
< |
GFaSeqGet* fastaSeqGet(GffObj& mrna) { |
324 |
< |
if (gseqpath.is_empty()) return NULL; |
325 |
< |
if (multiGSeq && mrna.gseq_id!=gseq_id) { |
326 |
< |
if (faseq!=NULL) { |
327 |
< |
delete faseq; |
328 |
< |
faseq=NULL; |
329 |
< |
} |
330 |
< |
faseq=openFastaSeq(mrna.getGSeqName()); |
331 |
< |
gseq_id=mrna.gseq_id; |
332 |
< |
} |
333 |
< |
return faseq; |
334 |
< |
} |
335 |
< |
|
336 |
< |
void adjust_stopcodon(GffObj& mrna, int adj, GList<GSeg>* seglst=NULL) { |
261 |
> |
int adjust_stopcodon(GffObj& gffrec, int adj, GList<GSeg>* seglst=NULL) { |
262 |
|
//adj>0 => extedn CDS, adj<0 => shrink CDS |
263 |
|
//when CDS is expanded, exons have to be checked too and |
264 |
|
// expanded accordingly if they had the same boundary |
265 |
< |
if (mrna.strand=='-') { |
266 |
< |
if ((int)mrna.CDstart>adj) { |
267 |
< |
mrna.CDstart-=adj; |
268 |
< |
if (adj>0 && mrna.exons.First()->start>=mrna.CDstart) { |
269 |
< |
mrna.exons.First()->start-=adj; |
270 |
< |
mrna.start=mrna.exons.First()->start; |
271 |
< |
mrna.covlen+=adj; |
265 |
> |
int realadj=0; |
266 |
> |
if (gffrec.strand=='-') { |
267 |
> |
if ((int)gffrec.CDstart>adj) { |
268 |
> |
|
269 |
> |
gffrec.CDstart-=adj; |
270 |
> |
realadj=adj; |
271 |
> |
if (adj<0) { //restore |
272 |
> |
if (gffrec.exons.First()->start==gffrec.CDstart+adj) { |
273 |
> |
gffrec.exons.First()->start-=adj; |
274 |
> |
gffrec.start=gffrec.exons.First()->start; |
275 |
> |
gffrec.covlen+=adj; |
276 |
> |
} |
277 |
> |
} |
278 |
> |
else if (gffrec.exons.First()->start>=gffrec.CDstart) { |
279 |
> |
gffrec.exons.First()->start-=adj; |
280 |
> |
gffrec.start=gffrec.exons.First()->start; |
281 |
> |
gffrec.covlen+=adj; |
282 |
|
} |
283 |
|
} |
284 |
|
} |
285 |
|
else { |
286 |
< |
mrna.CDend+=adj; |
287 |
< |
if (adj>0 && mrna.exons.Last()->end<=mrna.CDend) { |
288 |
< |
mrna.exons.Last()->end+=adj; |
289 |
< |
mrna.end=mrna.exons.Last()->end; |
290 |
< |
mrna.covlen+=adj; |
286 |
> |
realadj=adj; |
287 |
> |
gffrec.CDend+=adj; |
288 |
> |
if (adj<0) {//restore |
289 |
> |
if (gffrec.exons.Last()->end==gffrec.CDend-adj) { |
290 |
> |
gffrec.exons.Last()->end+=adj; |
291 |
> |
gffrec.end=gffrec.exons.Last()->end; |
292 |
> |
gffrec.covlen+=adj; |
293 |
> |
} |
294 |
> |
} |
295 |
> |
else if (gffrec.exons.Last()->end<=gffrec.CDend) { |
296 |
> |
gffrec.exons.Last()->end+=adj; |
297 |
> |
gffrec.end=gffrec.exons.Last()->end; |
298 |
> |
gffrec.covlen+=adj; |
299 |
|
} |
300 |
|
} |
301 |
|
if (seglst!=NULL) seglst->Last()->end+=adj; |
302 |
+ |
return realadj; |
303 |
|
} |
304 |
|
|
305 |
< |
void process_mRNA(GffObj& mrna) { |
306 |
< |
GStr gene(mrna.getGene()); |
307 |
< |
GStr defline(mrna.getID()); |
308 |
< |
if (!gene.is_empty() && strcmp(gene.chars(),mrna.getID())!=0) { |
309 |
< |
int* isonum=isoCounter.Find(gene.chars()); |
310 |
< |
if (isonum==NULL) { |
311 |
< |
isonum=new int(1); |
312 |
< |
isoCounter.Add(mrna.getGene(),isonum); |
305 |
> |
bool process_transcript(GFastaDb& gfasta, GffObj& gffrec) { |
306 |
> |
//returns true if the transcript passed the filter |
307 |
> |
char* gname=gffrec.getGeneName(); |
308 |
> |
if (gname==NULL) gname=gffrec.getGeneID(); |
309 |
> |
GStr defline(gffrec.getID()); |
310 |
> |
if (f_out && !fmtGTF) { |
311 |
> |
const char* tn=NULL; |
312 |
> |
if ((tn=gffrec.getAttr("transcript_name"))!=NULL) { |
313 |
> |
gffrec.addAttr("Name", tn); |
314 |
> |
gffrec.removeAttr("transcript_name"); |
315 |
> |
} |
316 |
> |
} |
317 |
> |
if (ensembl_convert && startsWith(gffrec.getID(), "ENS")) { |
318 |
> |
const char* tn=gffrec.getTrackName(); |
319 |
> |
gffrec.addAttr("type", tn); |
320 |
> |
//bool is_gene=false; |
321 |
> |
bool is_pseudo=false; |
322 |
> |
if (strcmp(tn, "protein_coding")==0 || gffrec.hasCDS()) |
323 |
> |
gffrec.setFeatureName("mRNA"); |
324 |
> |
else { |
325 |
> |
if (strcmp(tn, "processed_transcript")==0) |
326 |
> |
gffrec.setFeatureName("proc_RNA"); |
327 |
> |
else { |
328 |
> |
//is_gene=endsWith(tn, "gene"); |
329 |
> |
is_pseudo=strifind(tn, "pseudo"); |
330 |
> |
if (is_pseudo) { |
331 |
> |
gffrec.setFeatureName("pseudo_RNA"); |
332 |
> |
} |
333 |
> |
else if (endsWith(tn, "RNA")) { |
334 |
> |
gffrec.setFeatureName(tn); |
335 |
> |
} else gffrec.setFeatureName("misc_RNA"); |
336 |
> |
} |
337 |
> |
} |
338 |
|
} |
339 |
< |
else (*isonum)++; |
340 |
< |
defline.appendfmt(" gene=%s", gene.chars()); |
341 |
< |
} |
339 |
> |
if (gname && strcmp(gname, gffrec.getID())!=0) { |
340 |
> |
int* isonum=isoCounter.Find(gname); |
341 |
> |
if (isonum==NULL) { |
342 |
> |
isonum=new int(1); |
343 |
> |
isoCounter.Add(gname,isonum); |
344 |
> |
} |
345 |
> |
else (*isonum)++; |
346 |
> |
defline.appendfmt(" gene=%s", gname); |
347 |
> |
} |
348 |
|
int seqlen=0; |
349 |
+ |
|
350 |
|
const char* tlabel=tracklabel; |
351 |
< |
if (tlabel==NULL) tlabel=mrna.getTrackName(); |
351 |
> |
if (tlabel==NULL) tlabel=gffrec.getTrackName(); |
352 |
|
//defline.appendfmt(" track:%s",tlabel); |
353 |
|
char* cdsnt = NULL; |
354 |
|
char* cdsaa = NULL; |
355 |
|
int aalen=0; |
356 |
< |
for (int i=1;i<mrna.exons.Count();i++) { |
357 |
< |
int ilen=mrna.exons[i]->start-mrna.exons[i-1]->end-1; |
358 |
< |
if (ilen>2000000) |
356 |
> |
for (int i=1;i<gffrec.exons.Count();i++) { |
357 |
> |
int ilen=gffrec.exons[i]->start-gffrec.exons[i-1]->end-1; |
358 |
> |
if (ilen>4000000) |
359 |
|
GMessage("Warning: very large intron (%d) for transcript %s\n", |
360 |
< |
ilen, mrna.getID()); |
360 |
> |
ilen, gffrec.getID()); |
361 |
|
if (ilen>maxintron) { |
362 |
< |
return; |
362 |
> |
return false; |
363 |
|
} |
364 |
|
} |
365 |
|
GList<GSeg> seglst(false,true); |
366 |
< |
fastaSeqGet(mrna); |
367 |
< |
if (spliceCheck && mrna.exons.Count()>1) { |
366 |
> |
GFaSeqGet* faseq=fastaSeqGet(gfasta, gffrec); |
367 |
> |
if (spliceCheck && gffrec.exons.Count()>1) { |
368 |
|
//check introns for splice site consensi ( GT-AG, GC-AG or AT-AC ) |
369 |
|
if (faseq==NULL) GError("Error: no genomic sequence available!\n"); |
370 |
< |
int glen=mrna.end-mrna.start+1; |
371 |
< |
const char* gseq=faseq->subseq(mrna.start, glen); |
372 |
< |
bool revcompl=(mrna.strand=='-'); |
370 |
> |
int glen=gffrec.end-gffrec.start+1; |
371 |
> |
const char* gseq=faseq->subseq(gffrec.start, glen); |
372 |
> |
bool revcompl=(gffrec.strand=='-'); |
373 |
|
bool ssValid=true; |
374 |
< |
for (int e=1;e<mrna.exons.Count();e++) { |
375 |
< |
const char* intron=gseq+mrna.exons[e-1]->end+1-mrna.start; |
376 |
< |
//debug |
401 |
< |
//uint istart=mrna.exons[e-1]->end+1; |
402 |
< |
//uint iend=mrna.exons[e]->start-1; |
403 |
< |
int intronlen=mrna.exons[e]->start-mrna.exons[e-1]->end-1; |
374 |
> |
for (int e=1;e<gffrec.exons.Count();e++) { |
375 |
> |
const char* intron=gseq+gffrec.exons[e-1]->end+1-gffrec.start; |
376 |
> |
int intronlen=gffrec.exons[e]->start-gffrec.exons[e-1]->end-1; |
377 |
|
GSpliceSite acceptorSite(intron,intronlen,true, revcompl); |
378 |
|
GSpliceSite donorSite(intron,intronlen, false, revcompl); |
379 |
|
//GMessage("%c intron %d-%d : %s .. %s\n", |
380 |
< |
// mrna.strand, istart, iend, donorSite.nt, acceptorSite.nt); |
380 |
> |
// gffrec.strand, istart, iend, donorSite.nt, acceptorSite.nt); |
381 |
|
if (acceptorSite=="AG") { // GT-AG or GC-AG |
382 |
|
if (!donorSite.canonicalDonor()) { |
383 |
|
ssValid=false;break; |
391 |
|
//GFREE(gseq); |
392 |
|
if (!ssValid) { |
393 |
|
if (verbose) |
394 |
< |
GMessage("Invalid splice sites found for '%s'\n",mrna.getID()); |
395 |
< |
return; //don't print this one! |
394 |
> |
GMessage("Invalid splice sites found for '%s'\n",gffrec.getID()); |
395 |
> |
return false; //don't print this one! |
396 |
|
} |
397 |
|
} |
398 |
|
|
399 |
|
bool trprint=true; |
400 |
+ |
int stopCodonAdjust=0; |
401 |
|
int mCDphase=0; |
402 |
< |
if (mrna.CDphase=='1' || mrna.CDphase=='2') |
403 |
< |
mCDphase = mrna.CDphase-'0'; |
402 |
> |
bool hasStop=false; |
403 |
> |
if (gffrec.CDphase=='1' || gffrec.CDphase=='2') |
404 |
> |
mCDphase = gffrec.CDphase-'0'; |
405 |
|
if (f_y!=NULL || f_x!=NULL || validCDSonly) { |
406 |
|
if (faseq==NULL) GError("Error: no genomic sequence provided!\n"); |
407 |
|
//if (protmap && fullCDSonly) { |
408 |
< |
//if (protmap && (fullCDSonly || (mrna.qlen>0 && mrna.qend==mrna.qlen))) { |
408 |
> |
//if (protmap && (fullCDSonly || (gffrec.qlen>0 && gffrec.qend==gffrec.qlen))) { |
409 |
> |
|
410 |
|
if (validCDSonly) { //make sure the stop codon is always included |
411 |
< |
adjust_stopcodon(mrna,3); |
411 |
> |
//adjust_stopcodon(gffrec,3); |
412 |
> |
stopCodonAdjust=adjust_stopcodon(gffrec,3); |
413 |
|
} |
414 |
|
int strandNum=0; |
415 |
|
int phaseNum=0; |
416 |
|
CDS_CHECK: |
417 |
< |
cdsnt=mrna.getSpliced(faseq, true, &seqlen,NULL,NULL,&seglst); |
417 |
> |
cdsnt=gffrec.getSpliced(faseq, true, &seqlen,NULL,NULL,&seglst); |
418 |
|
if (cdsnt==NULL) trprint=false; |
419 |
|
if (validCDSonly) { |
420 |
|
cdsaa=translateDNA(cdsnt, aalen, seqlen); |
421 |
|
char* p=strchr(cdsaa,'.'); |
422 |
< |
bool hasStop=false; |
422 |
> |
hasStop=false; |
423 |
|
if (p!=NULL) { |
424 |
|
if (p-cdsaa>=aalen-2) { //stop found as the last codon |
425 |
|
*p='0';//remove it |
427 |
|
if (aalen-2==p-cdsaa) { |
428 |
|
//previous to last codon is the stop codon |
429 |
|
//so correct the CDS stop accordingly |
430 |
< |
adjust_stopcodon(mrna,-3, &seglst); |
430 |
> |
adjust_stopcodon(gffrec,-3, &seglst); |
431 |
> |
stopCodonAdjust=0; //clear artificial stop adjustment |
432 |
|
seqlen-=3; |
433 |
|
cdsnt[seqlen]=0; |
434 |
|
} |
442 |
|
//in-frame stop codon found |
443 |
|
if (altPhases && phaseNum<3) { |
444 |
|
phaseNum++; |
445 |
< |
mrna.CDphase = '0'+((mCDphase+phaseNum)%3); |
445 |
> |
gffrec.CDphase = '0'+((mCDphase+phaseNum)%3); |
446 |
|
GFREE(cdsaa); |
447 |
|
goto CDS_CHECK; |
448 |
|
} |
449 |
< |
if (mrna.exons.Count()==1 && bothStrands) { |
449 |
> |
if (gffrec.exons.Count()==1 && bothStrands) { |
450 |
|
strandNum++; |
451 |
|
phaseNum=0; |
452 |
|
if (strandNum<2) { |
453 |
|
GFREE(cdsaa); |
454 |
< |
mrna.strand = (mrna.strand=='-') ? '+':'-'; |
454 |
> |
gffrec.strand = (gffrec.strand=='-') ? '+':'-'; |
455 |
|
goto CDS_CHECK; //repeat the CDS check for a different frame |
456 |
|
} |
457 |
|
} |
458 |
< |
if (verbose) GMessage("In-frame STOP found for '%s'\n",mrna.getID()); |
458 |
> |
if (verbose) GMessage("In-frame STOP found for '%s'\n",gffrec.getID()); |
459 |
|
} //has in-frame STOP |
460 |
|
if (fullCDSonly) { |
461 |
|
if (!hasStop || cdsaa[0]!='M') trprint=false; |
465 |
|
if (!trprint) { |
466 |
|
GFREE(cdsnt); |
467 |
|
GFREE(cdsaa); |
468 |
< |
return; |
468 |
> |
return false; |
469 |
|
} |
470 |
< |
if (f_out!=NULL) { |
471 |
< |
if (fmtGTF) mrna.printGtf(f_out,tracklabel); |
472 |
< |
else mrna.printGff(f_out, tracklabel); |
473 |
< |
} |
470 |
> |
if (stopCodonAdjust>0 && !hasStop) { |
471 |
> |
//restore stop codon location |
472 |
> |
adjust_stopcodon(gffrec, -stopCodonAdjust, &seglst); |
473 |
> |
if (cdsnt!=NULL && seqlen>0) { |
474 |
> |
seqlen-=stopCodonAdjust; |
475 |
> |
cdsnt[seqlen]=0; |
476 |
> |
} |
477 |
> |
if (cdsaa!=NULL) aalen--; |
478 |
> |
} |
479 |
> |
|
480 |
|
if (f_y!=NULL) { //CDS translation fasta output requested |
481 |
|
//char* |
482 |
|
if (cdsaa==NULL) { //translate now if not done before |
483 |
|
cdsaa=translateDNA(cdsnt, aalen, seqlen); |
484 |
|
} |
485 |
< |
if (fullattr && mrna.attrs!=NULL) { |
485 |
> |
if (fullattr && gffrec.attrs!=NULL) { |
486 |
|
//append all attributes found for each transcripts |
487 |
< |
for (int i=0;i<mrna.attrs->Count();i++) { |
487 |
> |
for (int i=0;i<gffrec.attrs->Count();i++) { |
488 |
|
defline.append(" "); |
489 |
< |
defline.append(mrna.getAttrName(i)); |
489 |
> |
defline.append(gffrec.getAttrName(i)); |
490 |
|
defline.append("="); |
491 |
< |
defline.append(mrna.getAttrValue(i)); |
491 |
> |
defline.append(gffrec.getAttrValue(i)); |
492 |
|
} |
493 |
|
} |
494 |
|
printFasta(f_y, defline, cdsaa, aalen); |
496 |
|
if (f_x!=NULL) { //CDS only |
497 |
|
if (writeExonSegs) { |
498 |
|
defline.append(" loc:"); |
499 |
< |
defline.append(mrna.getGSeqName()); |
500 |
< |
defline.appendfmt("(%c)",mrna.strand); |
499 |
> |
defline.append(gffrec.getGSeqName()); |
500 |
> |
defline.appendfmt("(%c)",gffrec.strand); |
501 |
|
//warning: not CDS coordinates are written here, but the exon ones |
502 |
< |
defline+=(int)mrna.start; |
502 |
> |
defline+=(int)gffrec.start; |
503 |
|
defline+=(char)'-'; |
504 |
< |
defline+=(int)mrna.end; |
504 |
> |
defline+=(int)gffrec.end; |
505 |
|
// -- here these are CDS substring coordinates on the spliced sequence: |
506 |
|
defline.append(" segs:"); |
507 |
|
for (int i=0;i<seglst.Count();i++) { |
511 |
|
defline+=(int)seglst[i]->end; |
512 |
|
} |
513 |
|
} |
514 |
< |
if (fullattr && mrna.attrs!=NULL) { |
514 |
> |
if (fullattr && gffrec.attrs!=NULL) { |
515 |
|
//append all attributes found for each transcript |
516 |
< |
for (int i=0;i<mrna.attrs->Count();i++) { |
516 |
> |
for (int i=0;i<gffrec.attrs->Count();i++) { |
517 |
|
defline.append(" "); |
518 |
< |
defline.append(mrna.getAttrName(i)); |
518 |
> |
defline.append(gffrec.getAttrName(i)); |
519 |
|
defline.append("="); |
520 |
< |
defline.append(mrna.getAttrValue(i)); |
520 |
> |
defline.append(gffrec.getAttrValue(i)); |
521 |
|
} |
522 |
|
} |
523 |
|
printFasta(f_x, defline, cdsnt, seqlen); |
528 |
|
uint cds_start=0; |
529 |
|
uint cds_end=0; |
530 |
|
seglst.Clear(); |
531 |
< |
char* exont=mrna.getSpliced(faseq, false, &seqlen, &cds_start, &cds_end, &seglst); |
531 |
> |
char* exont=gffrec.getSpliced(faseq, false, &seqlen, &cds_start, &cds_end, &seglst); |
532 |
|
if (exont!=NULL) { |
533 |
< |
if (mrna.CDstart>0) { |
533 |
> |
if (gffrec.CDstart>0) { |
534 |
|
defline.appendfmt(" CDS=%d-%d", cds_start, cds_end); |
535 |
|
} |
536 |
|
if (writeExonSegs) { |
537 |
|
defline.append(" loc:"); |
538 |
< |
defline.append(mrna.getGSeqName()); |
538 |
> |
defline.append(gffrec.getGSeqName()); |
539 |
|
defline+=(char)'|'; |
540 |
< |
defline+=(int)mrna.start; |
540 |
> |
defline+=(int)gffrec.start; |
541 |
|
defline+=(char)'-'; |
542 |
< |
defline+=(int)mrna.end; |
542 |
> |
defline+=(int)gffrec.end; |
543 |
|
defline+=(char)'|'; |
544 |
< |
defline+=(char)mrna.strand; |
544 |
> |
defline+=(char)gffrec.strand; |
545 |
|
defline.append(" exons:"); |
546 |
< |
for (int i=0;i<mrna.exons.Count();i++) { |
546 |
> |
for (int i=0;i<gffrec.exons.Count();i++) { |
547 |
|
if (i>0) defline.append(","); |
548 |
< |
defline+=(int)mrna.exons[i]->start; |
548 |
> |
defline+=(int)gffrec.exons[i]->start; |
549 |
|
defline.append("-"); |
550 |
< |
defline+=(int)mrna.exons[i]->end; |
550 |
> |
defline+=(int)gffrec.exons[i]->end; |
551 |
|
} |
552 |
|
defline.append(" segs:"); |
553 |
|
for (int i=0;i<seglst.Count();i++) { |
557 |
|
defline+=(int)seglst[i]->end; |
558 |
|
} |
559 |
|
} |
560 |
< |
if (fullattr && mrna.attrs!=NULL) { |
560 |
> |
if (fullattr && gffrec.attrs!=NULL) { |
561 |
|
//append all attributes found for each transcripts |
562 |
< |
for (int i=0;i<mrna.attrs->Count();i++) { |
562 |
> |
for (int i=0;i<gffrec.attrs->Count();i++) { |
563 |
|
defline.append(" "); |
564 |
< |
defline.append(mrna.getAttrName(i)); |
564 |
> |
defline.append(gffrec.getAttrName(i)); |
565 |
|
defline.append("="); |
566 |
< |
defline.append(mrna.getAttrValue(i)); |
566 |
> |
defline.append(gffrec.getAttrValue(i)); |
567 |
|
} |
568 |
|
} |
569 |
|
printFasta(f_w, defline, exont, seqlen); |
570 |
|
GFREE(exont); |
571 |
|
} |
572 |
|
} //writing f_w (spliced exons) |
573 |
< |
return; |
573 |
> |
return true; |
574 |
|
} |
575 |
|
|
576 |
|
void openfw(FILE* &f, GArgs& args, char opt) { |
588 |
|
#define FWCLOSE(fh) if (fh!=NULL && fh!=stdout) fclose(fh) |
589 |
|
#define FRCLOSE(fh) if (fh!=NULL && fh!=stdin) fclose(fh) |
590 |
|
|
591 |
+ |
void printGff3Header(FILE* f, GArgs& args) { |
592 |
+ |
fprintf(f, "# "); |
593 |
+ |
args.printCmdLine(f); |
594 |
+ |
fprintf(f, "##gff-version 3\n"); |
595 |
+ |
//for (int i=0;i<gseqdata.Count();i++) { |
596 |
+ |
// |
597 |
+ |
//} |
598 |
+ |
} |
599 |
+ |
|
600 |
+ |
bool validateGffRec(GffObj* gffrec, GList<GffObj>* gfnew) { |
601 |
+ |
if (reftbl.Count()>0) { |
602 |
+ |
GStr refname(gffrec->getRefName()); |
603 |
+ |
RefTran* rt=reftbl.Find(refname.chars()); |
604 |
+ |
if (rt==NULL && refname[-2]=='.' && isdigit(refname[-1])) { |
605 |
+ |
//try removing the version |
606 |
+ |
refname.cut(-2); |
607 |
+ |
//GMessage("[DEBUG] Trying ref name '%s'...\n", refname.chars()); |
608 |
+ |
rt=reftbl.Find(refname.chars()); |
609 |
+ |
} |
610 |
+ |
if (rt) { |
611 |
+ |
gffrec->setRefName(rt->new_name); |
612 |
+ |
} |
613 |
+ |
else return false; //discard, ref seq not in the given translation table |
614 |
+ |
} |
615 |
+ |
if (mRNAOnly && gffrec->isDiscarded()) { |
616 |
+ |
//discard generic "gene" or "locus" features with no other detailed subfeatures |
617 |
+ |
//GMessage("Warning: discarding %s GFF generic gene/locus container %s\n",m->getID()); |
618 |
+ |
return false; |
619 |
+ |
} |
620 |
+ |
/* |
621 |
+ |
if (gffrec->exons.Count()==0 && gffrec->children.Count()==0)) { |
622 |
+ |
//a non-mRNA feature with no subfeatures |
623 |
+ |
//just so we get some sequence functions working, add a dummy "exon"-like subfeature here |
624 |
+ |
//--this could be a single "pseudogene" entry or another genomic region without exons |
625 |
+ |
// |
626 |
+ |
gffrec->addExon(gffrec->start,gffrec->end); |
627 |
+ |
} |
628 |
+ |
*/ |
629 |
+ |
if (rfltGSeq!=NULL) { //filter by gseqName |
630 |
+ |
if (strcmp(gffrec->getGSeqName(),rfltGSeq)!=0) { |
631 |
+ |
return false; |
632 |
+ |
} |
633 |
+ |
} |
634 |
+ |
if (rfltStrand>0 && gffrec->strand !=rfltStrand) { |
635 |
+ |
return false; |
636 |
+ |
} |
637 |
+ |
//check coordinates |
638 |
+ |
if (rfltStart!=0 || rfltEnd!=MAX_UINT) { |
639 |
+ |
if (rfltWithin) { |
640 |
+ |
if (gffrec->start<rfltStart || gffrec->end>rfltEnd) { |
641 |
+ |
return false; |
642 |
+ |
} |
643 |
+ |
} |
644 |
+ |
else { |
645 |
+ |
if (gffrec->start>rfltEnd || gffrec->end<rfltStart) { |
646 |
+ |
return false; |
647 |
+ |
} |
648 |
+ |
} |
649 |
+ |
} |
650 |
+ |
if (multiExon && gffrec->exons.Count()<=1) { |
651 |
+ |
return false; |
652 |
+ |
} |
653 |
+ |
if (wCDSonly && gffrec->CDstart==0) { |
654 |
+ |
return false; |
655 |
+ |
} |
656 |
+ |
if (ensembl_convert && startsWith(gffrec->getID(), "ENS")) { |
657 |
+ |
//keep track of chr|gene_id data -- coordinate range |
658 |
+ |
char* geneid=gffrec->getGeneID(); |
659 |
+ |
if (geneid!=NULL) { |
660 |
+ |
GeneInfo* ginfo=gene_ids.Find(geneid); |
661 |
+ |
if (ginfo==NULL) {//first time seeing this gene ID |
662 |
+ |
GeneInfo* geneinfo=new GeneInfo(gffrec, ensembl_convert); |
663 |
+ |
gene_ids.Add(geneid, geneinfo); |
664 |
+ |
if (gfnew!=NULL) gfnew->Add(geneinfo->gf); |
665 |
+ |
} |
666 |
+ |
else ginfo->update(gffrec); |
667 |
+ |
} |
668 |
+ |
} |
669 |
+ |
return true; |
670 |
+ |
} |
671 |
+ |
|
672 |
+ |
|
673 |
|
int main(int argc, char * const argv[]) { |
674 |
< |
GArgs args(argc, argv, "vOUNHWCVMNSXTDAPRZFGg:i:r:s:t:a:b:o:w:x:y:MINCOV=MINPID="); |
675 |
< |
int e; |
676 |
< |
if ((e=args.isError())>0) |
677 |
< |
GError("%s\nInvalid argument: %s\n", USAGE, argv[e]); |
678 |
< |
if (args.getOpt('h')!=NULL) GError("%s\n", USAGE); |
679 |
< |
debugMode=(args.getOpt('D')!=NULL); |
674 |
> |
GArgs args(argc, argv, |
675 |
> |
"debug;merge;cluster-only;help;MINCOV=MINPID=hvOUNHWCVJMKQNSXTDAPRZFGLEm:g:i:r:s:t:a:b:o:w:x:y:d:"); |
676 |
> |
args.printError(USAGE, true); |
677 |
> |
if (args.getOpt('h') || args.getOpt("help")) { |
678 |
> |
GMessage("%s",USAGE); |
679 |
> |
exit(1); |
680 |
> |
} |
681 |
> |
debugMode=(args.getOpt("debug")!=NULL); |
682 |
|
mRNAOnly=(args.getOpt('O')==NULL); |
683 |
< |
sortByLoc=(args.getOpt('S')!=NULL); |
683 |
> |
//sortByLoc=(args.getOpt('S')!=NULL); |
684 |
|
addDescr=(args.getOpt('A')!=NULL); |
685 |
|
verbose=(args.getOpt('v')!=NULL); |
686 |
|
wCDSonly=(args.getOpt('C')!=NULL); |
688 |
|
altPhases=(args.getOpt('H')!=NULL); |
689 |
|
fmtGTF=(args.getOpt('T')!=NULL); //switch output format to GTF |
690 |
|
bothStrands=(args.getOpt('B')!=NULL); |
691 |
< |
fullCDSonly=(args.getOpt('M')!=NULL); |
691 |
> |
fullCDSonly=(args.getOpt('J')!=NULL); |
692 |
|
spliceCheck=(args.getOpt('N')!=NULL); |
693 |
+ |
bool matchAllIntrons=(args.getOpt('K')==NULL); |
694 |
+ |
bool fuzzSpan=(args.getOpt('Q')!=NULL); |
695 |
+ |
if (args.getOpt('M') || args.getOpt("merge")) { |
696 |
+ |
doCluster=true; |
697 |
+ |
doCollapseRedundant=true; |
698 |
+ |
} |
699 |
+ |
else { |
700 |
+ |
if (!matchAllIntrons || fuzzSpan) { |
701 |
+ |
GMessage("%s",USAGE); |
702 |
+ |
GMessage("Error: -K or -Q options require -M/--merge option!\n"); |
703 |
+ |
exit(1); |
704 |
+ |
} |
705 |
+ |
} |
706 |
+ |
if (args.getOpt("cluster-only")) { |
707 |
+ |
doCluster=true; |
708 |
+ |
doCollapseRedundant=false; |
709 |
+ |
if (!matchAllIntrons || fuzzSpan) { |
710 |
+ |
GMessage("%s",USAGE); |
711 |
+ |
GMessage("Error: -K or -Q options have no effect with --cluster-only.\n"); |
712 |
+ |
exit(1); |
713 |
+ |
} |
714 |
+ |
} |
715 |
|
//protmap=(args.getOpt('P')!=NULL); |
716 |
|
if (fullCDSonly) validCDSonly=true; |
717 |
+ |
if (verbose) { |
718 |
+ |
fprintf(stderr, "Command line was:\n"); |
719 |
+ |
args.printCmdLine(stderr); |
720 |
+ |
} |
721 |
+ |
|
722 |
|
fullattr=(args.getOpt('F')!=NULL); |
723 |
|
if (args.getOpt('G')==NULL) |
724 |
|
noExonAttr=!fullattr; |
726 |
|
noExonAttr=true; |
727 |
|
fullattr=true; |
728 |
|
} |
729 |
+ |
ensembl_convert=(args.getOpt('L')!=NULL); |
730 |
+ |
if (ensembl_convert) { |
731 |
+ |
fullattr=true; |
732 |
+ |
noExonAttr=false; |
733 |
+ |
//sortByLoc=true; |
734 |
+ |
} |
735 |
+ |
|
736 |
|
mergeCloseExons=(args.getOpt('Z')!=NULL); |
737 |
|
multiExon=(args.getOpt('U')!=NULL); |
738 |
|
writeExonSegs=(args.getOpt('W')!=NULL); |
739 |
|
tracklabel=args.getOpt('t'); |
740 |
< |
gseqpath=args.getOpt('g'); |
741 |
< |
if (!gseqpath.is_empty()) { |
742 |
< |
sortByLoc=true; |
641 |
< |
int c=fileExists(gseqpath.chars()); |
642 |
< |
if (c==1) |
643 |
< |
multiGSeq=true;//directory given |
644 |
< |
else if (c==2) { //single file given |
645 |
< |
if (gseqpath.rindex(".cidx")==gseqpath.length()-5) { |
646 |
< |
//cdbyank index given |
647 |
< |
gcdb=new GCdbYank(gseqpath.chars()); |
648 |
< |
if (fileExists(gcdb->getDbName())) { |
649 |
< |
gcdbfa=gcdb->getDbName(); |
650 |
< |
} else { |
651 |
< |
gcdbfa=gseqpath; |
652 |
< |
gcdbfa.chomp(".cidx"); |
653 |
< |
if (!fileExists(gcdbfa.chars())) |
654 |
< |
GError("Error: cannot locate the fasta file for index %s !\n",gseqpath.chars()); |
655 |
< |
} |
656 |
< |
multiGSeq=true; |
657 |
< |
} |
658 |
< |
else { //must be just the one single fasta file |
659 |
< |
multiGSeq=false; |
660 |
< |
faseq=new GFaSeqGet(gseqpath.chars(), false); // WARNING: does not validate FASTA line-len uniformity! |
661 |
< |
} |
662 |
< |
}//fasta file given |
663 |
< |
else GError("Error: cannot find genomic sequence path %s !\n",gseqpath.chars()); |
664 |
< |
} // -g option |
740 |
> |
GFastaDb gfasta(args.getOpt('g')); |
741 |
> |
//if (gfasta.fastaPath!=NULL) |
742 |
> |
// sortByLoc=true; //enforce sorting by chromosome/contig |
743 |
|
GStr s=args.getOpt('i'); |
744 |
|
if (!s.is_empty()) maxintron=s.asInt(); |
745 |
+ |
|
746 |
+ |
FILE* f_repl=NULL; |
747 |
+ |
s=args.getOpt('d'); |
748 |
+ |
if (!s.is_empty()) { |
749 |
+ |
if (s=="-") f_repl=stdout; |
750 |
+ |
else { |
751 |
+ |
f_repl=fopen(s.chars(), "w"); |
752 |
+ |
if (f_repl==NULL) GError("Error creating file %s\n", s.chars()); |
753 |
+ |
} |
754 |
+ |
} |
755 |
+ |
|
756 |
|
rfltWithin=(args.getOpt('R')!=NULL); |
757 |
|
s=args.getOpt('r'); |
758 |
|
if (!s.is_empty()) { |
759 |
|
s.trim(); |
760 |
+ |
if (s[0]=='+' || s[0]=='-') { |
761 |
+ |
rfltStrand=s[0]; |
762 |
+ |
s.cut(0,1); |
763 |
+ |
} |
764 |
|
int isep=s.index(':'); |
765 |
< |
if (isep>=0) { //gseq name given |
766 |
< |
if (isep>0) rfltGSeq=Gstrdup((s.substr(0,isep)).chars()); |
765 |
> |
if (isep>0) { //gseq name given |
766 |
> |
if (rfltStrand==0 && (s[isep-1]=='+' || s[isep-1]=='-')) { |
767 |
> |
isep--; |
768 |
> |
rfltStrand=s[isep]; |
769 |
> |
s.cut(isep,1); |
770 |
> |
} |
771 |
> |
if (isep>0) |
772 |
> |
rfltGSeq=Gstrdup((s.substr(0,isep)).chars()); |
773 |
|
s.cut(0,isep+1); |
774 |
|
} |
775 |
|
GStr gsend; |
776 |
+ |
char slast=s[s.length()-1]; |
777 |
+ |
if (rfltStrand==0 && (slast=='+' || slast=='-')) { |
778 |
+ |
s.chomp(slast); |
779 |
+ |
rfltStrand=slast; |
780 |
+ |
} |
781 |
|
if (s.index("..")>=0) gsend=s.split(".."); |
782 |
< |
else gsend=s.split('-'); |
782 |
> |
else gsend=s.split('-'); |
783 |
|
if (!s.is_empty()) rfltStart=(uint)s.asInt(); |
784 |
< |
if (!gsend.is_empty()) rfltEnd=(uint)gsend.asInt(); |
784 |
> |
if (!gsend.is_empty()) { |
785 |
> |
rfltEnd=(uint)gsend.asInt(); |
786 |
> |
if (rfltEnd==0) rfltEnd=MAX_UINT; |
787 |
> |
} |
788 |
> |
|
789 |
|
} //gseq/range filtering |
790 |
|
else { |
791 |
|
if (rfltWithin) |
792 |
|
GError("Error: option -R doesn't make sense without -r!\n"); |
793 |
|
} |
794 |
+ |
s=args.getOpt('m'); |
795 |
+ |
if (!s.is_empty()) { |
796 |
+ |
FILE* ft=fopen(s,"r"); |
797 |
+ |
if (ft==NULL) GError("Error opening reference table: %s\n",s.chars()); |
798 |
+ |
loadRefTable(ft, reftbl); |
799 |
+ |
fclose(ft); |
800 |
+ |
} |
801 |
|
s=args.getOpt('s'); |
802 |
|
if (!s.is_empty()) { |
803 |
< |
FILE* fsize=fopen(s,"r"); |
804 |
< |
if (fsize==NULL) GError("Error opening info file: %s\n",s.chars()); |
805 |
< |
loadSeqInfo(fsize, seqinfo); |
806 |
< |
} |
807 |
< |
/* |
808 |
< |
openfw(fgtfok, args, 'a'); |
694 |
< |
openfw(fgtfbad, args, 'b'); |
695 |
< |
*/ |
803 |
> |
FILE* fsize=fopen(s,"r"); |
804 |
> |
if (fsize==NULL) GError("Error opening info file: %s\n",s.chars()); |
805 |
> |
loadSeqInfo(fsize, seqinfo); |
806 |
> |
fclose(fsize); |
807 |
> |
} |
808 |
> |
|
809 |
|
openfw(f_out, args, 'o'); |
810 |
|
//if (f_out==NULL) f_out=stdout; |
811 |
< |
if (gseqpath.is_empty() && (validCDSonly || spliceCheck || args.getOpt('w')!=NULL || args.getOpt('x')!=NULL || args.getOpt('y')!=NULL)) |
811 |
> |
if (gfasta.fastaPath==NULL && (validCDSonly || spliceCheck || args.getOpt('w')!=NULL || args.getOpt('x')!=NULL || args.getOpt('y')!=NULL)) |
812 |
|
GError("Error: -g option is required for options -w, -x, -y, -V, -N, -M !\n"); |
813 |
|
|
814 |
|
openfw(f_w, args, 'w'); |
816 |
|
openfw(f_y, args, 'y'); |
817 |
|
if (f_y!=NULL || f_x!=NULL) wCDSonly=true; |
818 |
|
//useBadCDS=useBadCDS || (fgtfok==NULL && fgtfbad==NULL && f_y==NULL && f_x==NULL); |
819 |
< |
GStr infile; |
820 |
< |
if (args.startNonOpt()) { |
821 |
< |
infile=args.nextNonOpt(); |
822 |
< |
//GMessage("Given file: %s\n",infile.chars()); |
823 |
< |
} |
824 |
< |
if (!infile.is_empty()) { |
825 |
< |
f_in=fopen(infile, "r"); |
826 |
< |
if (f_in==NULL) |
827 |
< |
GError("Cannot open input file %s!\n",infile.chars()); |
828 |
< |
} |
829 |
< |
else |
830 |
< |
f_in=stdin; |
831 |
< |
|
719 |
< |
//GffReader* gfreader=new GffReader(f_in,true); |
720 |
< |
GffReader* gfreader=new GffReader(f_in, mRNAOnly, sortByLoc); |
721 |
< |
gfreader->readAll(fullattr, mergeCloseExons, noExonAttr); |
722 |
< |
if (debugMode) GMessage("[D] gff data loaded, now processing each entry\n"); |
723 |
< |
//GList<GffObj> mrnas(true,true,false); |
724 |
< |
for (int m=0;m<gfreader->gflst.Count();m++) { |
725 |
< |
GffObj* mrna=gfreader->gflst[m]; |
726 |
< |
if (mrna->exons.Count()==0) { |
727 |
< |
//a non-mRNA feature with no subfeatures |
728 |
< |
//yet just so we get some sequence functions working, add the dummy "exon" here |
729 |
< |
mrna->addExon(mrna->start,mrna->end); |
730 |
< |
} |
731 |
< |
if (rfltGSeq!=NULL) { //filter by gseqName |
732 |
< |
if (strcmp(mrna->getGSeqName(),rfltGSeq)!=0) { |
733 |
< |
delete mrna; |
734 |
< |
gfreader->gflst.Forget(m); |
735 |
< |
continue; |
736 |
< |
} |
737 |
< |
} |
738 |
< |
//check coordinates |
739 |
< |
if (rfltStart!=0 || rfltEnd!=MAX_UINT) { |
740 |
< |
if (rfltWithin) { |
741 |
< |
if (mrna->start<rfltStart || mrna->end>rfltEnd) { |
742 |
< |
delete mrna; |
743 |
< |
gfreader->gflst.Forget(m); |
744 |
< |
continue; |
819 |
> |
|
820 |
> |
int numfiles = args.startNonOpt(); |
821 |
> |
//GList<GffObj> gfkept(false,true); //unsorted, free items on delete |
822 |
> |
int out_counter=0; //number of records printed |
823 |
> |
while (true) { |
824 |
> |
GStr infile; |
825 |
> |
if (numfiles) { |
826 |
> |
infile=args.nextNonOpt(); |
827 |
> |
if (infile.is_empty()) break; |
828 |
> |
if (infile=="-") { f_in=stdin; infile="stdin"; } |
829 |
> |
else |
830 |
> |
if ((f_in=fopen(infile, "r"))==NULL) |
831 |
> |
GError("Error: cannot open input file %s!\n",infile.chars()); |
832 |
|
} |
833 |
< |
} |
834 |
< |
else { |
835 |
< |
if (mrna->start>rfltEnd || mrna->end<rfltStart) { |
836 |
< |
delete mrna; |
837 |
< |
gfreader->gflst.Forget(m); |
838 |
< |
continue; |
833 |
> |
else |
834 |
> |
infile="-"; |
835 |
> |
GffLoader gffloader(infile.chars()); |
836 |
> |
gffloader.transcriptsOnly=mRNAOnly; |
837 |
> |
gffloader.fullAttributes=fullattr; |
838 |
> |
gffloader.noExonAttrs=noExonAttr; |
839 |
> |
gffloader.mergeCloseExons=mergeCloseExons; |
840 |
> |
gffloader.showWarnings=(args.getOpt('E')!=NULL); |
841 |
> |
gffloader.load(g_data, &validateGffRec, doCluster, doCollapseRedundant, |
842 |
> |
matchAllIntrons, fuzzSpan); |
843 |
> |
if (doCluster) |
844 |
> |
collectLocusData(g_data); |
845 |
> |
if (numfiles==0) break; |
846 |
> |
} |
847 |
> |
|
848 |
> |
GStr loctrack("gffcl"); |
849 |
> |
if (tracklabel) loctrack=tracklabel; |
850 |
> |
g_data.setSorted(&gseqCmpName); |
851 |
> |
if (doCluster) { |
852 |
> |
//grouped in loci |
853 |
> |
for (int g=0;g<g_data.Count();g++) { |
854 |
> |
GenomicSeqData* gdata=g_data[g]; |
855 |
> |
for (int l=0;l<gdata->loci.Count();l++) { |
856 |
> |
GffLocus& loc=*(gdata->loci[l]); |
857 |
> |
//check all non-replaced transcripts in this locus: |
858 |
> |
int numvalid=0; |
859 |
> |
int idxfirstvalid=-1; |
860 |
> |
for (int i=0;i<loc.rnas.Count();i++) { |
861 |
> |
GffObj& t=*(loc.rnas[i]); |
862 |
> |
GTData* tdata=(GTData*)(t.uptr); |
863 |
> |
if (tdata->replaced_by!=NULL) { |
864 |
> |
if (f_repl && (t.udata & 8)==0) { |
865 |
> |
//t.udata|=8; |
866 |
> |
fprintf(f_repl, "%s", t.getID()); |
867 |
> |
GTData* rby=tdata; |
868 |
> |
while (rby->replaced_by!=NULL) { |
869 |
> |
fprintf(f_repl," => %s", rby->replaced_by->getID()); |
870 |
> |
rby->rna->udata|=8; |
871 |
> |
rby=(GTData*)(rby->replaced_by->uptr); |
872 |
> |
} |
873 |
> |
fprintf(f_repl, "\n"); |
874 |
> |
} |
875 |
> |
continue; |
876 |
> |
} |
877 |
> |
if (process_transcript(gfasta, t)) { |
878 |
> |
t.udata|=4; //tag it as valid |
879 |
> |
numvalid++; |
880 |
> |
if (idxfirstvalid<0) idxfirstvalid=i; |
881 |
> |
} |
882 |
|
} |
883 |
+ |
|
884 |
+ |
if (f_out && numvalid>0) { |
885 |
+ |
GStr locname("RLOC_"); |
886 |
+ |
locname.appendfmt("%08d",loc.locus_num); |
887 |
+ |
if (!fmtGTF) { |
888 |
+ |
if (out_counter==0) |
889 |
+ |
printGff3Header(f_out, args); |
890 |
+ |
fprintf(f_out,"%s\t%s\tlocus\t%d\t%d\t.\t%c\t.\tID=%s;locus=%s", |
891 |
+ |
loc.rnas[0]->getGSeqName(), loctrack.chars(), loc.start, loc.end, loc.strand, |
892 |
+ |
locname.chars(), locname.chars()); |
893 |
+ |
//const char* loc_gname=loc.getGeneName(); |
894 |
+ |
if (loc.gene_names.Count()>0) { //print all gene names associated to this locus |
895 |
+ |
fprintf(f_out, ";genes=%s",loc.gene_names.First()->name.chars()); |
896 |
+ |
for (int i=1;i<loc.gene_names.Count();i++) { |
897 |
+ |
fprintf(f_out, ",%s",loc.gene_names[i]->name.chars()); |
898 |
+ |
} |
899 |
+ |
} |
900 |
+ |
if (loc.gene_ids.Count()>0) { //print all GeneIDs names associated to this locus |
901 |
+ |
fprintf(f_out, ";geneIDs=%s",loc.gene_ids.First()->name.chars()); |
902 |
+ |
for (int i=1;i<loc.gene_ids.Count();i++) { |
903 |
+ |
fprintf(f_out, ",%s",loc.gene_ids[i]->name.chars()); |
904 |
+ |
} |
905 |
+ |
} |
906 |
+ |
fprintf(f_out, ";transcripts=%s",loc.rnas[idxfirstvalid]->getID()); |
907 |
+ |
for (int i=idxfirstvalid+1;i<loc.rnas.Count();i++) { |
908 |
+ |
fprintf(f_out, ",%s",loc.rnas[i]->getID()); |
909 |
+ |
} |
910 |
+ |
fprintf(f_out, "\n"); |
911 |
+ |
} |
912 |
+ |
//now print all valid, non-replaced transcripts in this locus: |
913 |
+ |
for (int i=0;i<loc.rnas.Count();i++) { |
914 |
+ |
GffObj& t=*(loc.rnas[i]); |
915 |
+ |
GTData* tdata=(GTData*)(t.uptr); |
916 |
+ |
if (tdata->replaced_by!=NULL || ((t.udata & 4)==0)) continue; |
917 |
+ |
t.addAttr("locus", locname.chars()); |
918 |
+ |
out_counter++; |
919 |
+ |
if (fmtGTF) t.printGtf(f_out, tracklabel); |
920 |
+ |
else { |
921 |
+ |
//print the parent first, if any |
922 |
+ |
if (t.parent!=NULL && ((t.parent->udata & 4)==0)) { |
923 |
+ |
GTData* pdata=(GTData*)(t.parent->uptr); |
924 |
+ |
if (pdata->geneinfo!=NULL) |
925 |
+ |
pdata->geneinfo->finalize(); |
926 |
+ |
t.parent->addAttr("locus", locname.chars()); |
927 |
+ |
t.parent->printGff(f_out, tracklabel); |
928 |
+ |
t.parent->udata|=4; |
929 |
+ |
} |
930 |
+ |
t.printGff(f_out, tracklabel); |
931 |
+ |
} |
932 |
+ |
} |
933 |
+ |
} //have valid transcripts to print |
934 |
+ |
}//for each locus |
935 |
+ |
if (f_out && !mRNAOnly) { |
936 |
+ |
//final pass through the non-transcripts, in case any of them were not printed |
937 |
+ |
//TODO: order broken, these should be interspersed among the rnas in the correct order! |
938 |
+ |
for (int m=0;m<gdata->gfs.Count();m++) { |
939 |
+ |
GffObj& t=*(gdata->gfs[m]); |
940 |
+ |
if ((t.udata&4)==0) { //never printed |
941 |
+ |
t.udata|=4; |
942 |
+ |
if (fmtGTF) t.printGtf(f_out, tracklabel); |
943 |
+ |
else t.printGff(f_out, tracklabel); |
944 |
+ |
} |
945 |
+ |
} //for each non-transcript |
946 |
|
} |
947 |
< |
} |
948 |
< |
if (multiExon && mrna->exons.Count()==1) { |
949 |
< |
delete mrna; |
950 |
< |
gfreader->gflst.Forget(m); |
951 |
< |
continue; |
952 |
< |
} |
953 |
< |
if (wCDSonly && mrna->CDstart==0) { |
954 |
< |
delete mrna; |
955 |
< |
gfreader->gflst.Forget(m); |
956 |
< |
continue; |
957 |
< |
} |
958 |
< |
/* if (validCDSonly && mrna->hasErrors) { |
959 |
< |
delete mrna; |
960 |
< |
gfreader->gflst.Forget(m); |
961 |
< |
continue; |
947 |
> |
} //for each genomic sequence |
948 |
> |
} |
949 |
> |
else { |
950 |
> |
//not grouped into loci, print the rnas with their parents, if any |
951 |
> |
int numvalid=0; |
952 |
> |
for (int g=0;g<g_data.Count();g++) { |
953 |
> |
GenomicSeqData* gdata=g_data[g]; |
954 |
> |
for (int m=0;m<gdata->rnas.Count();m++) { |
955 |
> |
GffObj& t=*(gdata->rnas[m]); |
956 |
> |
GTData* tdata=(GTData*)(t.uptr); |
957 |
> |
if (tdata->replaced_by!=NULL) continue; |
958 |
> |
if (process_transcript(gfasta, t)) { |
959 |
> |
t.udata|=4; //tag it as valid |
960 |
> |
numvalid++; |
961 |
> |
if (f_out) { |
962 |
> |
if (tdata->geneinfo) tdata->geneinfo->finalize(); |
963 |
> |
out_counter++; |
964 |
> |
if (fmtGTF) t.printGtf(f_out, tracklabel); |
965 |
> |
else { |
966 |
> |
if (out_counter==1) |
967 |
> |
printGff3Header(f_out, args); |
968 |
> |
//print the parent first, if any |
969 |
> |
if (t.parent!=NULL && ((t.parent->udata & 4)==0)) { |
970 |
> |
GTData* pdata=(GTData*)(t.parent->uptr); |
971 |
> |
if (pdata->geneinfo!=NULL) |
972 |
> |
pdata->geneinfo->finalize(); |
973 |
> |
t.parent->printGff(f_out, tracklabel); |
974 |
> |
t.parent->udata|=4; |
975 |
> |
} |
976 |
> |
t.printGff(f_out, tracklabel); |
977 |
> |
} |
978 |
> |
}//GFF/GTF output requested |
979 |
> |
} //valid transcript |
980 |
> |
} //for each rna |
981 |
> |
if (f_out && !mRNAOnly) { |
982 |
> |
//final pass through the non-transcripts, in case any of them were not printed |
983 |
> |
//TODO: order broken, these should be interspersed among the rnas in the correct order! |
984 |
> |
for (int m=0;m<gdata->gfs.Count();m++) { |
985 |
> |
GffObj& t=*(gdata->gfs[m]); |
986 |
> |
if ((t.udata&4)==0) { //never printed |
987 |
> |
t.udata|=4; |
988 |
> |
if (fmtGTF) t.printGtf(f_out, tracklabel); |
989 |
> |
else t.printGff(f_out, tracklabel); |
990 |
> |
} |
991 |
> |
} //for each non-transcript |
992 |
|
} |
993 |
< |
*/ |
771 |
< |
process_mRNA(*mrna); |
772 |
< |
delete mrna; |
773 |
< |
gfreader->gflst.Forget(m); |
993 |
> |
} //for each genomic seq |
994 |
|
} |
995 |
< |
// M_END: |
776 |
< |
delete gfreader; |
995 |
> |
if (f_repl && f_repl!=stdout) fclose(f_repl); |
996 |
|
seqinfo.Clear(); |
997 |
< |
if (faseq!=NULL) delete faseq; |
998 |
< |
if (gcdb!=NULL) delete gcdb; |
997 |
> |
//if (faseq!=NULL) delete faseq; |
998 |
> |
//if (gcdb!=NULL) delete gcdb; |
999 |
|
GFREE(rfltGSeq); |
1000 |
|
FRCLOSE(f_in); |
1001 |
|
FWCLOSE(f_out); |