ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/scripts/gffmanip
Revision: 24
Committed: Tue Jul 26 21:46:39 2011 UTC (8 years, 2 months ago) by gpertea
File size: 18716 byte(s)
Log Message:
Line User Rev File contents
1 gpertea 23 #!/usr/bin/perl
2     use strict;
3     use Getopt::Std;
4     use FindBin;use lib $FindBin::Bin;
5    
6     my $usage = q/Usage:
7     gffmanip [-v] [-o outrecords.txt] [{-r|-c|-e} <rangedata>] [-C] [-t <track>]\
8     [-i [<attr>:]<ids.txt>] [-a <attr>=<value>[,...]] [-f <feature1>[,..]] \
9     [-k <attrlist>] [-l [{-T|-L|-l <attrlist>}] [<gff\/gtf_file>]
10    
11     Filtering options:
12     -C output only those transcripts having CDS features
13     -f output only the features matching any of the strings in the given
14     comma-delimited list (e.g. -f 'mRNA,exon,CDS')
15     -v show records that do NOT match the given filtering options (like grep's -v)
16     -r show only records overlapping any of the intervals given in <rangedata>,
17     which has the format:
18     [<strand>]<chr>[:<start>-<end>[,<start2>-<end2>,...]]
19     -c same as -r, but only shows records fully contained in the interval(s) given
20     -e outputs records whose exon\/CDS segments overlap the interval(s) given
21     -i only shows records whose ID match any of the entries in the file <ids.txt>;
22     unless another attribute is given as a prefix (<attr>:), the ID attribute
23     is used for GFF3 input, and transcript_id for GTF
24     -a only shows records having the <attr> attribute with the value <value>;
25     multiple attribute\/value pairs can be given (comma delimited) and a record
26     will be shown if there is at least one attribute match
27    
28     Output options:
29     -K keep all attributes ; by default only the core GFF\/GTF attributes are
30     shown (transcript_id, gene_id, gene_name, ID, Parent, Name, gene_name)
31     -k keep only the specified non-core attributes in the output records
32     -T output GTF format (default is GFF3)
33     -L output only the list of transcript IDs found in the file
34     (and matching the filtering options)
35     -l output the value of all attributes in <attrlist> for each record
36     (tab delimited if more than one attribute is given)
37     -t replace the 2nd column with the given <track> text
38     /;
39     umask 0002;
40     getopts('vTKACLf:r:l:c:e:i:a:t:k:o:') || die($usage."\n");
41     die("$usage\n") unless @ARGV>0;
42     die("Only one input file is expected (or '-' for stdin)!\n") unless @ARGV==1;
43     my $input_gff=$ARGV[0];
44     my $outfile=$Getopt::Std::opt_o;
45     my %ignoredFeatures;
46     @ignoredFeatures{qw(intron cds_start cds_stop start_codon stop_codon start stop cdsstart cdsstop)}=();
47     my %oattrs; #hash with attrs to print
48     @oattrs{qw(ID transcript_id Parent gene_id gene gene_name)}=();
49     my %attrflt; # attribute=>value filter
50     my %idflt; # hash with IDs to keep
51     my %atab; # hash with attribute names whose values will be listed as tab delimited
52     my @atab_cols;
53     if ($outfile) {
54     open(OUTF, '>'.$outfile) || die("Error creating output file $outfile\n");
55     select(OUTF);
56     }
57     # --
58     my ($fltinvert, $printGTF, $out_track)=
59     ($Getopt::Std::opt_v, $Getopt::Std::opt_T, $Getopt::Std::opt_t);
60     my $only_ifCDS=$Getopt::Std::opt_v;
61     my ($idfile, $avflt, $featlist, $oattrlist)=
62     ($Getopt::Std::opt_i, $Getopt::Std::opt_a, $Getopt::Std::opt_f, $Getopt::Std::opt_k);
63     my $keepAll=$Getopt::Std::opt_K || $Getopt::Std::opt_A || (lc($oattrlist) eq 'a' || lc($oattrlist) eq 'all');
64     my ($r_range, $c_range, $e_range)=
65     ($Getopt::Std::opt_r, $Getopt::Std::opt_c, $Getopt::Std::opt_e);
66     my $check_range;
67     if ($r_range) {
68     $check_range=$r_range;
69     die("Error: options -r, -c and -e are mutually exclusive!\n") if ($c_range || $e_range);
70     }
71     elsif ($c_range) {
72     $check_range=$c_range;
73     die("Error: options -r, -c and -e are mutually exclusive!\n") if ($r_range || $e_range);
74     }
75     elsif ($e_range) {
76     $check_range=$e_range;
77     die("Error: options -r, -c and -e are mutually exclusive!\n") if ($r_range || $c_range);
78     }
79     my ($flt_chr, $flt_strand);
80     my @flt_intv;
81     if ($check_range) {
82     ($flt_chr, my $ck_rlst)=split(/\:/,$check_range);
83     #die("$usage Incorrect format for the interval list!\n") unless $flt_chr && $ck_rlst;
84     my $flt_strand=substr($flt_chr,0,1);
85     if ($flt_strand eq '-' || $flt_strand eq '+') {
86     substr($flt_chr,0,1)='';
87     }
88     else {
89     $flt_strand=undef;
90     my $e=chop($flt_chr);
91     if ($e eq '-' || $e eq '+') {
92     $flt_strand=$e;
93     }
94     else { $flt_chr.=$e; }
95     } #no strand
96     my @flt_intv;
97     if ($ck_rlst) {
98     my @rdata=map { [split(/[\-\.]+/)] } (split(/[\,\;\s]+/,$ck_rlst));
99     foreach my $d (@rdata) {
100     ($$d[0], $$d[1])=($$d[1], $$d[0]) if $$d[0]>$$d[1];
101     }
102     @flt_intv = sort { $a->[0] <=> $b->[0] } @rdata;
103     }
104     }
105     #my $range=$flt_chr.':'.$ck_rex[0]->[0].'-'.$ex[-1]->[1];
106    
107     $oattrlist='' if $keepAll;
108     my $tab_attrs=$Getopt::Std::opt_l;
109     $tab_attrs='ID' if $Getopt::Std::opt_L;
110     if ($tab_attrs) {
111     @atab_cols=split(/\,/, $tab_attrs);
112     @atab{@atab_cols}=();
113     }
114     my $idfileattr; # attribute to use for ID list filtering instead of ID/transcript_id
115     if ($idfile) {
116     unless (-f $idfile) {
117     my ($a,$f)=(split(/\:/,$idfile));
118     ($idfile, $idfileattr)=($f,$a) if $f;
119     }
120     my $idf;
121     if ($idfile eq '-') {
122     open($idf, "<&=STDIN") || die("Error: couldn't alias STDIN. $!\n");
123     }
124     else {
125     open($idf, $idfile) || die("Error: cannot open $idfile $!\n");
126     }
127     while (<$idf>) {
128     my ($id)=(m/(\S+)/);
129     $idflt{$id}=1 if $id;
130     }
131     close($idf);
132     }
133     if ($avflt) {
134     my @avl=split(/\,/, $avflt);
135     foreach my $avpair (@avl) {
136     my ($attr, $value)=split(/\s*=\s*/,$avpair,2);
137     $value=~s/[" ]+$//;$value=~s/^[" ]+//;
138     $attrflt{$attr}=$value;
139     }
140     }
141     my %featflt; # feature list filter -- if $featlist was given
142     if ($featlist) {
143     #add to the list of attributes to be kept
144     my @fl=split(/\,/, $featlist);
145     @featflt{@fl}=(); # this is actually an input filter
146     }
147     if ($oattrlist) {
148     #add to the list of attributes to be kept
149     my @al=split(/\,/, $oattrlist);
150     @oattrs{@al}=();
151     }
152     # ----
153     my %gffrecs; # recID => [ chr, strand, feat_type, \%attrs, fstart, fend, [@exons], [@cds], isgff3, rejected, track, subfeat, fscore ]
154     # 0 1 2 3 4 5 6 7 8 9 10 , 11, 12
155     # recID has the prefix '<chr>|' which should be removed before output
156     my $gffh;
157     if ($input_gff eq '-') {
158     open($gffh, "<&=STDIN") || die("Error: couldn't alias STDIN $!\n");
159     }
160     else {
161     open($gffh, $input_gff) || die("Error opening file $input_gff $!\n");
162     }
163    
164     loadGff($gffh, \%gffrecs);
165     my @sorted_recs=sort sortByLoc keys(%gffrecs);
166     processGffRecs(\%gffrecs, \@sorted_recs);
167    
168     # --
169     if ($outfile) {
170     select(STDOUT);
171     close(OUTF);
172     }
173    
174     #************ Subroutines **************
175     sub sortByLoc {
176     my $da=$gffrecs{$a};
177     my $db=$gffrecs{$b};
178     if ($$da[0] eq $$db[0]) {
179     return ($$da[4]==$$db[4]) ? $$da[5] <=> $$db[5] : $$da[4] <=> $$db[4] ;
180     }
181     else { return $$da[0] cmp $$db[0] ; }
182     }
183    
184     sub checkOvlExons {
185     my ($a, $b, $rx)=@_;
186     return 0 if ($a>$$rx[-1]->[1] || $b<$$rx[0]->[0]); # not overlapping the whole exon chain
187     foreach my $x (@$rx) {
188     return 1 if ($a<=$$x[1] && $b>=$$x[0]);
189     return 0 if $b<$$x[0];
190     }
191     }
192    
193     sub checkWithinExons {
194     my ($a, $b, $rx)=@_; #checks if interval $a-$b is contained in any @$rx interval
195     return 0 if ($a>$$rx[-1]->[1] || $b<$$rx[0]->[0]); # not overlapping the whole exon chain
196     foreach my $x (@$rx) {
197     return 1 if ($a>=$$x[0] && $b<=$$x[1]);
198     return 0 if $b<$$x[0];
199     }
200     }
201    
202    
203     sub loadGff {
204     my ($gffhandle, $recs)=@_;
205     while (<$gffhandle>) {
206     next if m/^\s*#/;
207     chomp;
208     my ($chr, $track, $f, $fstart, $fend, $fscore, $strand, $frame, $lnum)=split(/\t/);
209     next unless $fstart>1 && $lnum;
210     $track=$out_track if $out_track;
211     next if exists($ignoredFeatures{lc($f)});
212     $f='exon' if $f=~m/exon/i;
213     $f='CDS' if $f=~m/^cds$/i;
214     my $rejected;
215     $rejected=1 if ($featlist && not exists($featflt{$f}));
216     #next if $f eq 'gene' || $f eq 'locus'; # Warning: skipping any 'gene' or 'locus' features, unconditionally
217     my $gff3_ID;
218     my $gff3_Parent;
219     my ($gname,$tdescr);
220     my %attrs;
221     ($fstart, $fend)=($fend, $fstart) if $fend<$fstart;
222     #$lnum=~s/"([^"]+)\;([^"]+)"/"$1.$2"/g; #protect ; within text between quotes
223     my @av=split(/\s*\;\s*/,$lnum);
224     ($gff3_ID)=($lnum=~m/\bID=([^;]+)/);
225     ($gff3_Parent)=($lnum=~m/\bParent=([^;]+)/);
226     my $isGFF3 = ($gff3_ID || $gff3_Parent);
227     if ($isGFF3) { # GFF format
228     $gff3_ID=~tr/"//d; #"
229     $gff3_Parent=~tr/"//d; #"
230     $gff3_Parent='' if ($f =~m/RNA/); # we really don't care about parent for RNA features
231     if ($gff3_ID && !$gff3_Parent) { #top level feature
232     foreach my $a (@av) {
233     my ($attr, $value)=split(/\s*=\s*/,$a,2);
234     $value=~s/[" ]+$//;$value=~s/^[" ]+//;
235     $attrs{$attr}=$value;
236     }
237     if ($f=~m/RNA/i || $f=~/gene/) {
238     # try to parse the description, if any
239    
240     if ($lnum=~m/\b(?:descr|tophit|info|product)\s*=\s*"?([^;"]+)/i) {
241     $tdescr=$1;
242     }
243     elsif ($lnum=~m/Name\s*=\s*"?([^;"]+)/) {
244     $tdescr=$1;
245     }
246     if ($lnum=~m/\bgene_name[\s=]+"?([^;"]+)/i) {
247     $gname=$1;
248     }
249     elsif ($lnum=~m/Name\s*=\s*"?([^;"]+)/) {
250     $gname=$1;
251     }
252     $tdescr='' if ($tdescr eq $gname);
253     $gname='' if $gname eq $gff3_ID;
254     }
255     die("Error: duplicate feature $gff3_ID on $chr\n") if (exists($recs->{"$chr|$gff3_ID"}));
256     my $recID="$chr|$gff3_ID";
257     $recs->{$recID} = [$chr, $strand, $f, {%attrs}, $fstart, $fend, [], [], $isGFF3, $rejected, $track, '', $fscore];
258     next;
259     } # parent/top-level feature
260     } #GFF
261     else { #GTF format
262     if ($f eq 'transcript') { # GTF with parent 'transcript' feature
263     foreach my $a (@av) {
264     my ($attr, $value)=split(/\s+"/,$a,2); #"
265     $value=~s/[" ]+$//;
266     if ($attr eq 'transcript_id') {
267     $attr='ID';
268     }
269     elsif ($attr eq 'gene_id') {
270     $attr='Name';
271     }
272     $attrs{$attr}=$value;
273     }
274     my $recID=$attrs{'ID'};
275     die("Error: cannot find transcript_id for GTF 'transcript' line:\n$_\n") unless $recID;
276     die("Error: duplicate feature $recID on $chr\n") if (exists($recs->{"$chr|$recID"}));
277     $recID=$chr.'|'.$recID;
278     $recs->{$recID} = [$chr, $strand, $f, {%attrs}, $fstart, $fend, [], [], $isGFF3, $rejected, $track, '', $fscore ];
279     next;
280     } # parent 'transcript' feature in GTF
281     }
282     # -------------- exon/CDS line here:
283     next if ($featlist && !exists($featflt{$f}));
284     my $recID;
285     ($gname, $tdescr)=();
286     if ($isGFF3) {
287     $recID=$gff3_Parent;
288     }
289     elsif ($lnum=~m/transcript_id[= ]+(['"\:\w\.\|\-]+)/) {
290     $recID=$1;
291     $recID=~tr/"//d; #"
292     }
293     elsif ($track=~/^jigsaw/ && $lnum=~m/^\d+$/) {
294     $recID=$chr.'.jsm.'.$lnum;
295     $gff3_Parent=$recID;
296     $isGFF3=1;
297     $f='CDS';
298     }
299     else {
300     die("Error: cannot parse locus/transcript name from input line:\n$_\n");
301     }
302     if (!$gname && $lnum=~m/gene_id[= ]+(['"\:\w\.\|\-]+)/) {
303     $gname=$1;
304     $gname=~tr/"//d; #"
305     }
306     $tdescr='' if index($recID, $tdescr)>=0;
307     $gname='' if index($recID, $gname)>=0;
308     $recID=$chr.'|'.$recID;
309     my $ld = $recs->{$recID};
310     if ($ld) { #existing entry
311     my $i=($f eq 'CDS') ? 7 : 6;
312     my ($lstart, $lend)=($$ld[4], $$ld[5]);
313     $$ld[4]=$fstart if $fstart<$lstart;
314     $$ld[5]=$fend if $fend>$lend;
315     push(@{$$ld[$i]}, [$fstart, $fend, $fscore, $frame]);
316     if ($f ne 'CDS') {
317     if ($$ld[11] && $$ld[11] ne $f) {
318     die("Error: multiple non-CDS subfeatures found for $recID ($$ld[11], $f)\n");
319     }
320     $$ld[11]=$f;
321     }
322     }
323     else { # first time seeing this locus/gene
324     # get the attributes from this first exon line
325     if ($gff3_Parent) { #has GFF3 Parent
326     foreach my $a (@av) {
327     my ($attr, $value)=split(/\s*=\s*/,$a,2); #"
328     $attr='ID' if $attr eq 'Parent';
329     $value=~s/[" ]+$//;$value=~s/^[" ]+//;
330     }
331     }
332     else { # GTF
333     foreach my $a (@av) {
334     my ($attr, $value)=split(/\s+"/,$a,2); #"
335     $value=~s/[" ]+$//;
336     if ($attr eq 'transcript_id') {
337     $attr='ID';
338     }
339     elsif ($attr eq 'gene_id') {
340     $attr='Name';
341     }
342     $attrs{$attr}=$value unless $attr=~m/^exon/;;
343     }
344     }
345     $recs->{$recID} = ($f eq 'CDS') ?
346     [$chr, $strand, $f, {%attrs}, $fstart, $fend, [], [[$fstart, $fend, $fscore, $frame]], $isGFF3, $rejected, $track, '','.' ] :
347     [$chr, $strand, $f, {%attrs}, $fstart, $fend, [[$fstart, $fend, $fscore, $frame]], [], $isGFF3, $rejected, $track, $f,'.' ] ;
348     # 0 1 2 3 4 5 6(exons) 7 (CDS) 8 9 10
349     }
350     }
351     close($gffh);
352     }
353    
354     sub processGffRecs {
355     #return if keys(%recs)==0;
356     my ($recs, $rlist)=@_;
357     my @recs_keys;
358     unless ($rlist) {
359     @recs_keys=keys(%$recs);
360     $rlist=\@recs_keys;
361     }
362     foreach my $recid (@$rlist) {
363     my $td=$$recs{$recid};
364     # 0 1 2 3 4 5 6 7 8 9 10 11 12
365     my ($chr, $strand, $feature, $attrs, $fstart, $fend, $er, $cr, $isGFF3, $rej, $track, $subf, $fscore) = @$td;
366     next if ($rej && !$fltinvert);
367     next if ($fltinvert && $featlist && !$rej);
368     # my ($mstart,$mend)=($fstart, $fend);
369     my $CDSonly=0; # set to 1 if only CDS segments were given
370     my $hasCDS=0;
371     my @ex;
372     my @cds;
373     #some records might lack exons, but have only CDS segments (e.g. mitochondrial genes)
374     if (@$er<1 && @$cr>0) {
375     @ex = sort { $a->[0] <=> $b->[0] } @$cr;
376     @cds=@ex;
377     $CDSonly=1;
378     $hasCDS=1;
379     }
380     else {
381     @ex = sort { $a->[0] <=> $b->[0] } @$er;
382     if (@$cr>0) { # sort cds segments too
383     @cds= sort { $a->[0] <=> $b->[0] } @$cr;
384     $hasCDS=1;
385     }
386     }
387     # --------------
388     # get the more accurate version of the start-end coords for the feature
389     my $covlen=0;
390     # map { $covlen+=$_->[1]-$_->[0]+1 } @ex;
391     my ($mstart, $mend) = ($ex[0]->[0], $ex[-1]->[1]);
392     my $gffid=$recid;
393     substr($gffid, 0, length($chr)+1)='';
394     die("Error: gffid ($gffid) not matching attrs\{ID\}(".$attrs-{'ID'}.")!\n")
395     unless ($gffid eq $attrs->{'ID'});
396     my $doprint=1;
397     #check all the filters
398     $doprint=0 if ($only_ifCDS && $hasCDS==0);
399     if ($doprint && $idfile) {
400     if ($idfileattr) {
401     $doprint=0 unless exists($idflt{$attrs->{$idfileattr}});
402     }
403     else {
404     $doprint=0 unless exists($idflt{$gffid});
405     }
406     }
407     if ($doprint && $avflt) {
408     my $avfound=0;
409     foreach my $a (keys(%attrflt)) {
410     if ($attrflt{$a} eq $attrs->{$a}) {
411     $avfound=1;
412     last;
413     }
414     }
415     $doprint=0 unless $avfound;
416     }
417    
418    
419     if ($doprint && $flt_chr && $flt_chr ne $chr) {
420     $doprint=0;
421     }
422     if ($doprint) {
423     if ($flt_strand && $flt_strand ne $strand) {
424     $doprint=0;
425     }
426     }
427     if ($doprint && @flt_intv>0) {
428     if ($r_range) { #t span overlap any range intervals
429     $doprint=1 if checkOvlExons($fstart, $fend, \@flt_intv);
430     }
431     elsif ($c_range) { #t span contained in any of the intervals
432     $doprint=1 if checkWithinExons($fstart, $fend, \@flt_intv);
433     }
434     elsif ($e_range) { # each exon checked for overlap of any interval
435     for my $ed (@ex) {
436     if (checkOvlExons($$ed[0], $$ed[1], \@flt_intv)) {
437     $doprint=1;
438     last;
439     }
440     }
441     }
442     }
443    
444     $doprint = ! $doprint if $fltinvert;
445     next unless $doprint;
446     #filter passed, print output
447     if ($tab_attrs) {
448     my @od;
449     foreach my $a (@atab_cols) {
450     if (uc($a) eq 'ID') {
451     push(@od, $gffid);
452     next;
453     }
454     push(@od, $attrs->{$a});
455     }
456     print join("\t",@od)."\n";
457     next;
458     }
459     #the core attributes ID, Name must be there already
460     unless (exists($attrs->{Name})) {
461     foreach my $name (qw(gene_name gene geneID geneId geneid locus loc ID)) {
462     if (exists($attrs->{$name})) {
463     $attrs->{Name}=$attrs->{$name};
464     last;
465     }
466     }
467     }
468     my ($tid, $tname)= (delete($attrs->{ID}), delete($attrs->{Name}));
469     my ($gene_name, $gene, $locus);
470     unless ($keepAll) {
471     if (exists($attrs->{gene_name})) {
472     $gene_name=delete($attrs->{gene_name});
473     $gene_name=undef if ($gene_name eq $tname && !exists($oattrs{gene_name}));
474     }
475     if (exists($attrs->{gene})) {
476     $gene=delete($attrs->{gene});
477     $gene=undef if ($gene eq $tname || $gene eq $gene_name) && !exists($oattrs{gene});
478     }
479     if (exists($attrs->{locus})) {
480     $locus=delete($attrs->{locus});
481     $locus=undef if ($locus eq $gene_name || $locus eq $gene_name) && !exists($oattrs{locus});
482     }
483     }
484     my $tattrs;
485     if ($printGTF) { #print GTF
486     $tattrs='transcript_id "'.$tid.'"; gene_id "'.$tname.'";';
487     $tattrs.=' gene_name "'.$gene_name.'";' if $gene_name;
488     $tattrs.=' gene "'.$gene.'";' if $gene;
489     $tattrs.=' locus "'.$locus.'";' if $locus;
490     if ($keepAll) {
491     foreach my $attr (keys(%$attrs)) {
492     $tattrs.=' '.$attr.' "'.$attrs->{$attr}.'";';
493     }
494     }
495     else { # only a subset of attributes will be shown
496     foreach my $attr (keys(%$attrs)) {
497     next unless exists($oattrs{$attr});
498     $tattrs.=' '.$attr.' "'.$attrs->{$attr}.'";';
499     }
500     }
501     }
502     else { #print GFF3
503     $tattrs='Parent='.$tid;
504     my $pattrs="ID=$tid;Name=$tname";
505     $pattrs.='gene_name='.$gene_name.';' if $gene_name;
506     $pattrs.='gene='.$gene.';' if $gene;
507     $pattrs.='locus='.$locus.';' if $locus;
508     if ($keepAll) {
509     foreach my $attr (keys(%$attrs)) {
510     my $val=$attrs->{$attr};
511     $pattrs.="$attr=$val;";
512     }
513     }
514     else { # only a subset of attributes will be shown
515     foreach my $attr (keys(%$attrs)) {
516     next unless exists($oattrs{$attr});
517     my $val=$attrs->{$attr};
518     $pattrs.="$attr=$val;";
519     }
520     }
521     print join("\t",$chr, $track, $feature, $fstart, $fend, $fscore, $strand, '.', $pattrs)."\n";
522     }
523     if ($CDSonly==0) { #write exons only when found in the input
524     foreach my $ed (@ex) {
525     print join("\t",$chr, $track, $subf, $$ed[0], $$ed[1], $$ed[2], $strand, $$ed[3], $tattrs)."\n";
526     }
527     }
528     foreach my $cd (@cds) {
529     print join("\t",$chr, $track, 'CDS', $$cd[0], $$cd[1], $$cd[2], $strand, $$cd[3], $tattrs)."\n";
530     }
531     } #for each stored transcript
532     }

Properties

Name Value
svn:executable *