1 |
gpertea |
17 |
#include "gff_utils.h" |
2 |
gpertea |
6 |
#include "GArgs.h" |
3 |
|
|
#include <ctype.h> |
4 |
gpertea |
17 |
// don't care about cdb compression |
5 |
|
|
//#ifdef ENABLE_COMPRESSION |
6 |
|
|
//#undef ENABLE_COMPRESSION |
7 |
|
|
//#endif |
8 |
|
|
//#include "GCdbYank.h" |
9 |
gpertea |
6 |
|
10 |
|
|
#define USAGE "Usage:\n\ |
11 |
gpertea |
17 |
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 |
gpertea |
6 |
\n\ |
18 |
|
|
Options:\n\ |
19 |
gpertea |
17 |
-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 |
gpertea |
6 |
for each of the mapped sequences:\n\ |
24 |
|
|
<seq-name> <seq-length> <seq-description>\n\ |
25 |
gpertea |
17 |
(useful for -A option with mRNA/EST/protein mappings)\n\ |
26 |
gpertea |
6 |
-i discard transcripts having an intron larger than <maxintron>\n\ |
27 |
|
|
-r only show transcripts crossing coordinate range <start>..<end>\n\ |
28 |
gpertea |
17 |
(on chromosome/contig <chr>, strand <strand> if provided)\n\ |
29 |
gpertea |
6 |
-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 |
gpertea |
17 |
-C coding only: discard mRNAs that have no CDS feature\n\ |
33 |
|
|
-F full GFF attribute preservation (all attributes are shown)\n\ |
34 |
gpertea |
6 |
-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 |
gpertea |
17 |
-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 |
gpertea |
6 |
\n\ |
39 |
gpertea |
17 |
-O process also non-transcript GFF records (by default non-transcript\n\ |
40 |
|
|
records are ignored)\n\ |
41 |
gpertea |
6 |
-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 |
gpertea |
17 |
-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 |
gpertea |
6 |
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 |
gpertea |
17 |
-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 |
gpertea |
6 |
coordinates projected onto the spliced sequence\n\ |
69 |
gpertea |
17 |
-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 |
gpertea |
6 |
|
82 |
gpertea |
17 |
|
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 |
gpertea |
6 |
FILE* ffasta=NULL; |
114 |
|
|
FILE* f_in=NULL; |
115 |
|
|
FILE* f_out=NULL; |
116 |
|
|
FILE* f_w=NULL; //fasta with spliced exons (transcripts) |
117 |
|
|
FILE* f_x=NULL; //fasta with spliced CDS |
118 |
|
|
FILE* f_y=NULL; //fasta with translated CDS |
119 |
|
|
bool wCDSonly=false; |
120 |
|
|
|
121 |
|
|
bool validCDSonly=false; // translation with no in-frame STOP |
122 |
|
|
bool bothStrands=false; //for single-exon mRNA validation, check the other strand too |
123 |
|
|
bool altPhases=false; //if original phase fails translation validation, |
124 |
|
|
//try the other 2 phases until one makes it |
125 |
|
|
bool mRNAOnly=true; |
126 |
|
|
bool spliceCheck=false; //only known splice-sites |
127 |
|
|
|
128 |
|
|
bool fullCDSonly=false; // starts with START, ends with STOP codon |
129 |
|
|
bool fullattr=false; |
130 |
gpertea |
17 |
//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 |
gpertea |
6 |
bool fmtGTF=false; |
141 |
|
|
bool addDescr=false; |
142 |
|
|
//bool protmap=false; |
143 |
|
|
bool multiExon=false; |
144 |
|
|
bool writeExonSegs=false; |
145 |
|
|
char* tracklabel=NULL; |
146 |
|
|
int maxintron=999000000; |
147 |
|
|
bool mergeCloseExons=false; |
148 |
|
|
//range filter: |
149 |
|
|
char* rfltGSeq=NULL; |
150 |
gpertea |
17 |
char rfltStrand=0; |
151 |
gpertea |
6 |
uint rfltStart=0; |
152 |
|
|
uint rfltEnd=MAX_UINT; |
153 |
|
|
bool rfltWithin=false; //check for full containment within given range |
154 |
|
|
bool noExonAttr=false; |
155 |
|
|
|
156 |
gpertea |
17 |
bool doCluster=false; |
157 |
|
|
bool doCollapseRedundant=false; |
158 |
gpertea |
6 |
|
159 |
gpertea |
17 |
GList<GenomicSeqData> g_data(true,true,true); //list of GFF records by genomic seq |
160 |
gpertea |
6 |
|
161 |
|
|
//hash with sequence info |
162 |
|
|
GHash<SeqInfo> seqinfo; |
163 |
|
|
GHash<int> isoCounter; //counts the valid isoforms |
164 |
gpertea |
17 |
GHash<RefTran> reftbl; |
165 |
|
|
GHash<GeneInfo> gene_ids; |
166 |
|
|
//min-max gene span associated to chr|gene_id (mostly for Ensembl conversion) |
167 |
gpertea |
6 |
|
168 |
|
|
bool debugMode=false; |
169 |
|
|
bool verbose=false; |
170 |
|
|
|
171 |
gpertea |
17 |
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 |
gpertea |
6 |
char* getSeqDescr(char* seqid) { |
217 |
|
|
static char charbuf[128]; |
218 |
|
|
if (seqinfo.Count()==0) return NULL; |
219 |
|
|
char* suf=rstrchr(seqid, '.'); |
220 |
|
|
if (suf!=NULL) *suf=0; |
221 |
|
|
SeqInfo* seqd=seqinfo.Find(seqid); |
222 |
|
|
if (suf!=NULL) *suf='.'; |
223 |
|
|
if (seqd!=NULL) { |
224 |
|
|
GStr s(seqd->descr); |
225 |
|
|
//cleanup some Uniref gunk |
226 |
|
|
if (s[0]=='[') { |
227 |
|
|
int r=s.index(']'); |
228 |
|
|
if (r>=0 && r<8 && isdigit(s[1])) |
229 |
|
|
s.remove(0,r+1); |
230 |
|
|
} |
231 |
|
|
if (s.length()>80) { |
232 |
|
|
int r=s.index(';'); |
233 |
|
|
if (r>5) s.cut(r); |
234 |
|
|
} |
235 |
|
|
if (s.length()>127) { |
236 |
|
|
s.cut(127); |
237 |
|
|
int r=s.rindex(' '); |
238 |
|
|
if (r>0) s.cut(r); |
239 |
|
|
} |
240 |
|
|
strcpy(charbuf, s.chars()); |
241 |
|
|
return charbuf; |
242 |
|
|
} |
243 |
|
|
else return NULL; |
244 |
|
|
} |
245 |
|
|
|
246 |
|
|
char* getSeqName(char* seqid) { |
247 |
|
|
static char charbuf[128]; |
248 |
|
|
char* suf=rstrchr(seqid, '.'); |
249 |
|
|
if (suf!=NULL) *suf=0; |
250 |
|
|
strcpy(charbuf, seqid); |
251 |
|
|
if (suf!=NULL) *suf='.'; |
252 |
|
|
return charbuf; |
253 |
|
|
} |
254 |
|
|
|
255 |
gpertea |
17 |
GFaSeqGet* fastaSeqGet(GFastaDb& gfasta, GffObj& gffrec) { |
256 |
|
|
if (gfasta.fastaPath==NULL) return NULL; |
257 |
|
|
return gfasta.fetch(gffrec.gseq_id); |
258 |
gpertea |
6 |
} |
259 |
|
|
|
260 |
|
|
|
261 |
gpertea |
17 |
int adjust_stopcodon(GffObj& gffrec, int adj, GList<GSeg>* seglst=NULL) { |
262 |
gpertea |
6 |
//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 |
gpertea |
17 |
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 |
gpertea |
6 |
} |
277 |
gpertea |
17 |
} |
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 |
gpertea |
6 |
} |
284 |
|
|
} |
285 |
|
|
else { |
286 |
gpertea |
17 |
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 |
gpertea |
6 |
} |
300 |
|
|
} |
301 |
|
|
if (seglst!=NULL) seglst->Last()->end+=adj; |
302 |
gpertea |
17 |
return realadj; |
303 |
gpertea |
6 |
} |
304 |
|
|
|
305 |
gpertea |
17 |
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 |
gpertea |
6 |
} |
339 |
gpertea |
17 |
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 |
gpertea |
6 |
int seqlen=0; |
349 |
gpertea |
17 |
|
350 |
gpertea |
6 |
const char* tlabel=tracklabel; |
351 |
gpertea |
17 |
if (tlabel==NULL) tlabel=gffrec.getTrackName(); |
352 |
gpertea |
6 |
//defline.appendfmt(" track:%s",tlabel); |
353 |
|
|
char* cdsnt = NULL; |
354 |
|
|
char* cdsaa = NULL; |
355 |
|
|
int aalen=0; |
356 |
gpertea |
17 |
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 |
gpertea |
6 |
GMessage("Warning: very large intron (%d) for transcript %s\n", |
360 |
gpertea |
17 |
ilen, gffrec.getID()); |
361 |
gpertea |
6 |
if (ilen>maxintron) { |
362 |
gpertea |
17 |
return false; |
363 |
gpertea |
6 |
} |
364 |
|
|
} |
365 |
|
|
GList<GSeg> seglst(false,true); |
366 |
gpertea |
17 |
GFaSeqGet* faseq=fastaSeqGet(gfasta, gffrec); |
367 |
|
|
if (spliceCheck && gffrec.exons.Count()>1) { |
368 |
gpertea |
6 |
//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 |
gpertea |
17 |
int glen=gffrec.end-gffrec.start+1; |
371 |
|
|
const char* gseq=faseq->subseq(gffrec.start, glen); |
372 |
|
|
bool revcompl=(gffrec.strand=='-'); |
373 |
gpertea |
6 |
bool ssValid=true; |
374 |
gpertea |
17 |
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 |
gpertea |
6 |
GSpliceSite acceptorSite(intron,intronlen,true, revcompl); |
378 |
|
|
GSpliceSite donorSite(intron,intronlen, false, revcompl); |
379 |
|
|
//GMessage("%c intron %d-%d : %s .. %s\n", |
380 |
gpertea |
17 |
// gffrec.strand, istart, iend, donorSite.nt, acceptorSite.nt); |
381 |
gpertea |
6 |
if (acceptorSite=="AG") { // GT-AG or GC-AG |
382 |
|
|
if (!donorSite.canonicalDonor()) { |
383 |
|
|
ssValid=false;break; |
384 |
|
|
} |
385 |
|
|
} |
386 |
|
|
else if (acceptorSite=="AC") { // |
387 |
|
|
if (donorSite!="AT") { ssValid=false; break; } |
388 |
|
|
} |
389 |
|
|
else { ssValid=false; break; } |
390 |
|
|
} |
391 |
|
|
//GFREE(gseq); |
392 |
|
|
if (!ssValid) { |
393 |
|
|
if (verbose) |
394 |
gpertea |
17 |
GMessage("Invalid splice sites found for '%s'\n",gffrec.getID()); |
395 |
|
|
return false; //don't print this one! |
396 |
gpertea |
6 |
} |
397 |
|
|
} |
398 |
|
|
|
399 |
|
|
bool trprint=true; |
400 |
gpertea |
17 |
int stopCodonAdjust=0; |
401 |
gpertea |
6 |
int mCDphase=0; |
402 |
gpertea |
17 |
bool hasStop=false; |
403 |
|
|
if (gffrec.CDphase=='1' || gffrec.CDphase=='2') |
404 |
|
|
mCDphase = gffrec.CDphase-'0'; |
405 |
gpertea |
6 |
if (f_y!=NULL || f_x!=NULL || validCDSonly) { |
406 |
|
|
if (faseq==NULL) GError("Error: no genomic sequence provided!\n"); |
407 |
|
|
//if (protmap && fullCDSonly) { |
408 |
gpertea |
17 |
//if (protmap && (fullCDSonly || (gffrec.qlen>0 && gffrec.qend==gffrec.qlen))) { |
409 |
|
|
|
410 |
gpertea |
6 |
if (validCDSonly) { //make sure the stop codon is always included |
411 |
gpertea |
17 |
//adjust_stopcodon(gffrec,3); |
412 |
|
|
stopCodonAdjust=adjust_stopcodon(gffrec,3); |
413 |
gpertea |
6 |
} |
414 |
|
|
int strandNum=0; |
415 |
|
|
int phaseNum=0; |
416 |
|
|
CDS_CHECK: |
417 |
gpertea |
17 |
cdsnt=gffrec.getSpliced(faseq, true, &seqlen,NULL,NULL,&seglst); |
418 |
gpertea |
6 |
if (cdsnt==NULL) trprint=false; |
419 |
|
|
if (validCDSonly) { |
420 |
|
|
cdsaa=translateDNA(cdsnt, aalen, seqlen); |
421 |
|
|
char* p=strchr(cdsaa,'.'); |
422 |
gpertea |
17 |
hasStop=false; |
423 |
gpertea |
6 |
if (p!=NULL) { |
424 |
|
|
if (p-cdsaa>=aalen-2) { //stop found as the last codon |
425 |
|
|
*p='0';//remove it |
426 |
|
|
hasStop=true; |
427 |
|
|
if (aalen-2==p-cdsaa) { |
428 |
|
|
//previous to last codon is the stop codon |
429 |
|
|
//so correct the CDS stop accordingly |
430 |
gpertea |
17 |
adjust_stopcodon(gffrec,-3, &seglst); |
431 |
|
|
stopCodonAdjust=0; //clear artificial stop adjustment |
432 |
gpertea |
6 |
seqlen-=3; |
433 |
|
|
cdsnt[seqlen]=0; |
434 |
|
|
} |
435 |
|
|
aalen=p-cdsaa; |
436 |
|
|
} |
437 |
|
|
else {//stop found before the last codon |
438 |
|
|
trprint=false; |
439 |
|
|
} |
440 |
|
|
}//stop codon found |
441 |
|
|
if (trprint==false) { //failed CDS validity check |
442 |
|
|
//in-frame stop codon found |
443 |
|
|
if (altPhases && phaseNum<3) { |
444 |
|
|
phaseNum++; |
445 |
gpertea |
17 |
gffrec.CDphase = '0'+((mCDphase+phaseNum)%3); |
446 |
gpertea |
6 |
GFREE(cdsaa); |
447 |
|
|
goto CDS_CHECK; |
448 |
|
|
} |
449 |
gpertea |
17 |
if (gffrec.exons.Count()==1 && bothStrands) { |
450 |
gpertea |
6 |
strandNum++; |
451 |
|
|
phaseNum=0; |
452 |
|
|
if (strandNum<2) { |
453 |
|
|
GFREE(cdsaa); |
454 |
gpertea |
17 |
gffrec.strand = (gffrec.strand=='-') ? '+':'-'; |
455 |
gpertea |
6 |
goto CDS_CHECK; //repeat the CDS check for a different frame |
456 |
|
|
} |
457 |
|
|
} |
458 |
gpertea |
17 |
if (verbose) GMessage("In-frame STOP found for '%s'\n",gffrec.getID()); |
459 |
gpertea |
6 |
} //has in-frame STOP |
460 |
|
|
if (fullCDSonly) { |
461 |
|
|
if (!hasStop || cdsaa[0]!='M') trprint=false; |
462 |
|
|
} |
463 |
|
|
} // CDS check requested |
464 |
|
|
} //translation or codon check/output was requested |
465 |
|
|
if (!trprint) { |
466 |
|
|
GFREE(cdsnt); |
467 |
|
|
GFREE(cdsaa); |
468 |
gpertea |
17 |
return false; |
469 |
gpertea |
6 |
} |
470 |
gpertea |
17 |
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 |
gpertea |
6 |
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 |
gpertea |
17 |
if (fullattr && gffrec.attrs!=NULL) { |
486 |
gpertea |
6 |
//append all attributes found for each transcripts |
487 |
gpertea |
17 |
for (int i=0;i<gffrec.attrs->Count();i++) { |
488 |
gpertea |
6 |
defline.append(" "); |
489 |
gpertea |
17 |
defline.append(gffrec.getAttrName(i)); |
490 |
gpertea |
6 |
defline.append("="); |
491 |
gpertea |
17 |
defline.append(gffrec.getAttrValue(i)); |
492 |
gpertea |
6 |
} |
493 |
|
|
} |
494 |
|
|
printFasta(f_y, defline, cdsaa, aalen); |
495 |
|
|
} |
496 |
|
|
if (f_x!=NULL) { //CDS only |
497 |
|
|
if (writeExonSegs) { |
498 |
|
|
defline.append(" loc:"); |
499 |
gpertea |
17 |
defline.append(gffrec.getGSeqName()); |
500 |
|
|
defline.appendfmt("(%c)",gffrec.strand); |
501 |
gpertea |
6 |
//warning: not CDS coordinates are written here, but the exon ones |
502 |
gpertea |
17 |
defline+=(int)gffrec.start; |
503 |
gpertea |
6 |
defline+=(char)'-'; |
504 |
gpertea |
17 |
defline+=(int)gffrec.end; |
505 |
gpertea |
6 |
// -- here these are CDS substring coordinates on the spliced sequence: |
506 |
|
|
defline.append(" segs:"); |
507 |
|
|
for (int i=0;i<seglst.Count();i++) { |
508 |
|
|
if (i>0) defline.append(","); |
509 |
|
|
defline+=(int)seglst[i]->start; |
510 |
|
|
defline.append("-"); |
511 |
|
|
defline+=(int)seglst[i]->end; |
512 |
|
|
} |
513 |
|
|
} |
514 |
gpertea |
17 |
if (fullattr && gffrec.attrs!=NULL) { |
515 |
gpertea |
6 |
//append all attributes found for each transcript |
516 |
gpertea |
17 |
for (int i=0;i<gffrec.attrs->Count();i++) { |
517 |
gpertea |
6 |
defline.append(" "); |
518 |
gpertea |
17 |
defline.append(gffrec.getAttrName(i)); |
519 |
gpertea |
6 |
defline.append("="); |
520 |
gpertea |
17 |
defline.append(gffrec.getAttrValue(i)); |
521 |
gpertea |
6 |
} |
522 |
|
|
} |
523 |
|
|
printFasta(f_x, defline, cdsnt, seqlen); |
524 |
|
|
} |
525 |
|
|
GFREE(cdsnt); |
526 |
|
|
GFREE(cdsaa); |
527 |
|
|
if (f_w!=NULL) { //write spliced exons |
528 |
|
|
uint cds_start=0; |
529 |
|
|
uint cds_end=0; |
530 |
|
|
seglst.Clear(); |
531 |
gpertea |
17 |
char* exont=gffrec.getSpliced(faseq, false, &seqlen, &cds_start, &cds_end, &seglst); |
532 |
gpertea |
6 |
if (exont!=NULL) { |
533 |
gpertea |
17 |
if (gffrec.CDstart>0) { |
534 |
gpertea |
6 |
defline.appendfmt(" CDS=%d-%d", cds_start, cds_end); |
535 |
|
|
} |
536 |
|
|
if (writeExonSegs) { |
537 |
|
|
defline.append(" loc:"); |
538 |
gpertea |
17 |
defline.append(gffrec.getGSeqName()); |
539 |
gpertea |
6 |
defline+=(char)'|'; |
540 |
gpertea |
17 |
defline+=(int)gffrec.start; |
541 |
gpertea |
6 |
defline+=(char)'-'; |
542 |
gpertea |
17 |
defline+=(int)gffrec.end; |
543 |
gpertea |
6 |
defline+=(char)'|'; |
544 |
gpertea |
17 |
defline+=(char)gffrec.strand; |
545 |
gpertea |
6 |
defline.append(" exons:"); |
546 |
gpertea |
17 |
for (int i=0;i<gffrec.exons.Count();i++) { |
547 |
gpertea |
6 |
if (i>0) defline.append(","); |
548 |
gpertea |
17 |
defline+=(int)gffrec.exons[i]->start; |
549 |
gpertea |
6 |
defline.append("-"); |
550 |
gpertea |
17 |
defline+=(int)gffrec.exons[i]->end; |
551 |
gpertea |
6 |
} |
552 |
|
|
defline.append(" segs:"); |
553 |
|
|
for (int i=0;i<seglst.Count();i++) { |
554 |
|
|
if (i>0) defline.append(","); |
555 |
|
|
defline+=(int)seglst[i]->start; |
556 |
|
|
defline.append("-"); |
557 |
|
|
defline+=(int)seglst[i]->end; |
558 |
|
|
} |
559 |
|
|
} |
560 |
gpertea |
17 |
if (fullattr && gffrec.attrs!=NULL) { |
561 |
gpertea |
6 |
//append all attributes found for each transcripts |
562 |
gpertea |
17 |
for (int i=0;i<gffrec.attrs->Count();i++) { |
563 |
gpertea |
6 |
defline.append(" "); |
564 |
gpertea |
17 |
defline.append(gffrec.getAttrName(i)); |
565 |
gpertea |
6 |
defline.append("="); |
566 |
gpertea |
17 |
defline.append(gffrec.getAttrValue(i)); |
567 |
gpertea |
6 |
} |
568 |
|
|
} |
569 |
|
|
printFasta(f_w, defline, exont, seqlen); |
570 |
|
|
GFREE(exont); |
571 |
|
|
} |
572 |
|
|
} //writing f_w (spliced exons) |
573 |
gpertea |
17 |
return true; |
574 |
gpertea |
6 |
} |
575 |
|
|
|
576 |
|
|
void openfw(FILE* &f, GArgs& args, char opt) { |
577 |
|
|
GStr s=args.getOpt(opt); |
578 |
|
|
if (!s.is_empty()) { |
579 |
|
|
if (s=='-') |
580 |
|
|
f=stdout; |
581 |
|
|
else { |
582 |
|
|
f=fopen(s,"w"); |
583 |
|
|
if (f==NULL) GError("Error creating file: %s\n", s.chars()); |
584 |
|
|
} |
585 |
|
|
} |
586 |
|
|
} |
587 |
|
|
|
588 |
|
|
#define FWCLOSE(fh) if (fh!=NULL && fh!=stdout) fclose(fh) |
589 |
|
|
#define FRCLOSE(fh) if (fh!=NULL && fh!=stdin) fclose(fh) |
590 |
|
|
|
591 |
gpertea |
17 |
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 |
gpertea |
6 |
int main(int argc, char * const argv[]) { |
674 |
gpertea |
17 |
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 |
gpertea |
6 |
mRNAOnly=(args.getOpt('O')==NULL); |
683 |
gpertea |
17 |
//sortByLoc=(args.getOpt('S')!=NULL); |
684 |
gpertea |
6 |
addDescr=(args.getOpt('A')!=NULL); |
685 |
|
|
verbose=(args.getOpt('v')!=NULL); |
686 |
|
|
wCDSonly=(args.getOpt('C')!=NULL); |
687 |
|
|
validCDSonly=(args.getOpt('V')!=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 |
gpertea |
17 |
fullCDSonly=(args.getOpt('J')!=NULL); |
692 |
gpertea |
6 |
spliceCheck=(args.getOpt('N')!=NULL); |
693 |
gpertea |
17 |
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 |
gpertea |
6 |
//protmap=(args.getOpt('P')!=NULL); |
716 |
|
|
if (fullCDSonly) validCDSonly=true; |
717 |
gpertea |
17 |
if (verbose) { |
718 |
|
|
fprintf(stderr, "Command line was:\n"); |
719 |
|
|
args.printCmdLine(stderr); |
720 |
|
|
} |
721 |
|
|
|
722 |
gpertea |
6 |
fullattr=(args.getOpt('F')!=NULL); |
723 |
|
|
if (args.getOpt('G')==NULL) |
724 |
|
|
noExonAttr=!fullattr; |
725 |
|
|
else { |
726 |
|
|
noExonAttr=true; |
727 |
|
|
fullattr=true; |
728 |
|
|
} |
729 |
gpertea |
17 |
ensembl_convert=(args.getOpt('L')!=NULL); |
730 |
|
|
if (ensembl_convert) { |
731 |
|
|
fullattr=true; |
732 |
|
|
noExonAttr=false; |
733 |
|
|
//sortByLoc=true; |
734 |
|
|
} |
735 |
|
|
|
736 |
gpertea |
6 |
mergeCloseExons=(args.getOpt('Z')!=NULL); |
737 |
|
|
multiExon=(args.getOpt('U')!=NULL); |
738 |
|
|
writeExonSegs=(args.getOpt('W')!=NULL); |
739 |
|
|
tracklabel=args.getOpt('t'); |
740 |
gpertea |
17 |
GFastaDb gfasta(args.getOpt('g')); |
741 |
|
|
//if (gfasta.fastaPath!=NULL) |
742 |
|
|
// sortByLoc=true; //enforce sorting by chromosome/contig |
743 |
gpertea |
6 |
GStr s=args.getOpt('i'); |
744 |
|
|
if (!s.is_empty()) maxintron=s.asInt(); |
745 |
gpertea |
17 |
|
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 |
gpertea |
6 |
rfltWithin=(args.getOpt('R')!=NULL); |
757 |
|
|
s=args.getOpt('r'); |
758 |
|
|
if (!s.is_empty()) { |
759 |
|
|
s.trim(); |
760 |
gpertea |
17 |
if (s[0]=='+' || s[0]=='-') { |
761 |
|
|
rfltStrand=s[0]; |
762 |
|
|
s.cut(0,1); |
763 |
|
|
} |
764 |
gpertea |
6 |
int isep=s.index(':'); |
765 |
gpertea |
17 |
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 |
gpertea |
6 |
s.cut(0,isep+1); |
774 |
|
|
} |
775 |
|
|
GStr gsend; |
776 |
gpertea |
17 |
char slast=s[s.length()-1]; |
777 |
|
|
if (rfltStrand==0 && (slast=='+' || slast=='-')) { |
778 |
|
|
s.chomp(slast); |
779 |
|
|
rfltStrand=slast; |
780 |
|
|
} |
781 |
gpertea |
6 |
if (s.index("..")>=0) gsend=s.split(".."); |
782 |
gpertea |
17 |
else gsend=s.split('-'); |
783 |
gpertea |
6 |
if (!s.is_empty()) rfltStart=(uint)s.asInt(); |
784 |
gpertea |
17 |
if (!gsend.is_empty()) { |
785 |
|
|
rfltEnd=(uint)gsend.asInt(); |
786 |
|
|
if (rfltEnd==0) rfltEnd=MAX_UINT; |
787 |
|
|
} |
788 |
|
|
|
789 |
gpertea |
6 |
} //gseq/range filtering |
790 |
|
|
else { |
791 |
|
|
if (rfltWithin) |
792 |
|
|
GError("Error: option -R doesn't make sense without -r!\n"); |
793 |
|
|
} |
794 |
gpertea |
17 |
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 |
gpertea |
6 |
s=args.getOpt('s'); |
802 |
|
|
if (!s.is_empty()) { |
803 |
gpertea |
17 |
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 |
gpertea |
6 |
openfw(f_out, args, 'o'); |
810 |
|
|
//if (f_out==NULL) f_out=stdout; |
811 |
gpertea |
17 |
if (gfasta.fastaPath==NULL && (validCDSonly || spliceCheck || args.getOpt('w')!=NULL || args.getOpt('x')!=NULL || args.getOpt('y')!=NULL)) |
812 |
gpertea |
6 |
GError("Error: -g option is required for options -w, -x, -y, -V, -N, -M !\n"); |
813 |
|
|
|
814 |
|
|
openfw(f_w, args, 'w'); |
815 |
|
|
openfw(f_x, args, 'x'); |
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 |
gpertea |
17 |
|
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 |
gpertea |
6 |
} |
833 |
gpertea |
17 |
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 |
gpertea |
6 |
} |
883 |
gpertea |
17 |
|
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 |
gpertea |
6 |
} |
947 |
gpertea |
17 |
} //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 |
gpertea |
6 |
} |
993 |
gpertea |
17 |
} //for each genomic seq |
994 |
gpertea |
6 |
} |
995 |
gpertea |
17 |
if (f_repl && f_repl!=stdout) fclose(f_repl); |
996 |
gpertea |
6 |
seqinfo.Clear(); |
997 |
gpertea |
17 |
//if (faseq!=NULL) delete faseq; |
998 |
|
|
//if (gcdb!=NULL) delete gcdb; |
999 |
gpertea |
6 |
GFREE(rfltGSeq); |
1000 |
|
|
FRCLOSE(f_in); |
1001 |
|
|
FWCLOSE(f_out); |
1002 |
|
|
FWCLOSE(f_w); |
1003 |
|
|
FWCLOSE(f_x); |
1004 |
|
|
FWCLOSE(f_y); |
1005 |
|
|
} |