| 1 |
#!/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 |
} |
