1 |
#!/usr/bin/perl |
2 |
use strict; |
3 |
my $usage=q/ |
4 |
gff_add_descr.pl <cdbyank_file.cidx> <gfflines..> |
5 |
|
6 |
Adds an extra "Descr" attribute to the main gff lines |
7 |
(mRNA or gene) of the input, by querying the ID |
8 |
from the <cdbyank_file.cidx> file. |
9 |
|
10 |
/; |
11 |
my %cache; # ID => descr |
12 |
my $cdbidx=shift(@ARGV) || die("$usage\nError: no cidx file given\n"); |
13 |
my $skipped=0; |
14 |
my $skip; |
15 |
while (<>) { |
16 |
if (m/^\s*#/) { print $_;next; } |
17 |
my @t=split('\t'); |
18 |
next unless($t[8]); #can't be a valid gff3 line |
19 |
my $descr; |
20 |
my $f=lc($t[2]); |
21 |
if ($f eq 'mrna' || $f eq 'gene') { |
22 |
if ($t[8]=~m/(?:Descr|Info)\s*=\s*"?([^;"]+)/i) { |
23 |
$descr=$1; |
24 |
if (length($descr)>6) { |
25 |
$skipped++; |
26 |
print $_; |
27 |
next; |
28 |
} |
29 |
$t[8]=~s/(?:Descr|Info)\s*=\s*"?[^;"]+//i; |
30 |
$t[8]=~tr/;/;/s; |
31 |
}#has descr/info already |
32 |
# -- retrieve the description here.. |
33 |
my ($id)=($t[8]=~m/\bID\s*=\s*"?([^;" ]+)/i); |
34 |
die("Error: no ID found for gff $f line: $_\n") unless $id; |
35 |
$id=~s/\.(\d+)$/.m$1/; |
36 |
$t[8]=~s/\bID\s*=\s*"?([^;" ]+)/ID=$id/; |
37 |
$id=~s/\.[a-z]*\d+$//; |
38 |
$descr=fetchDescr($id); |
39 |
chomp($t[8]); |
40 |
$t[8].=';descr="'.$descr.'"'; |
41 |
$_=join("\t",@t)."\n"; |
42 |
} #mRNA/gene line |
43 |
else { # |
44 |
s/(Parent=[\w\|]+)\.(\d+)$/$1.m$2/; |
45 |
} |
46 |
print $_; |
47 |
} #while |
48 |
|
49 |
sub fetchDescr { |
50 |
my $id=$_[0]; |
51 |
my $def=$cache{$id}; |
52 |
return $def if $def; |
53 |
#print STDERR "fetching: cdbyank -a '$id' -F $cdbidx\n"; |
54 |
$def=`cdbyank -a '$id' -F $cdbidx`; |
55 |
chomp($def); |
56 |
$def=~s/^\S+\s*//; #remove first token (the ID) |
57 |
return $def; |
58 |
} |