ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/scripts/gff_rangeflt.pl
Revision: 71
Committed: Fri Sep 23 19:23:49 2011 UTC (7 years, 11 months ago) by gpertea
File size: 1815 byte(s)
Log Message:
gridx fixed to export PERL5LIB and PYTHONPATH

Line File contents
1 #!/usr/bin/perl
2 use strict;
3 use Getopt::Std;
4 getopts('T') || die($usage."\n");
5
6 my $usage=q/
7 Usage:
8 gff_rangeflt.pl [-T] <chr>[<strand>]:<start_coord>-<end_coord> <input_gff..>
9
10 Use the -T option to not check specifically for exon overlap, but any overlap
11 of the transcript range (including intron-only overlaps)
12 /;
13
14 my $r=shift(@ARGV) || die("$usage\n");
15 die("$usage\n") if $r=~m/^\-/;
16 my ($chr, $cstart, $cend)=($r=~m/^([^:]+):(\d+)[-\.]+(\d+)/);
17 ($cstart, $cend)=($cend,$cstart) if $cstart>$cend;
18 die("$usage\n") unless $cstart>1 && $cend>$cstart;
19 my $cstrand=chop($chr);
20 unless ($cstrand eq '-' || $cstrand eq '+') {
21 $chr.=$cstrand; #restore, it wasn't a strand there
22 $cstrand='';
23 }
24 my @gffbuf; #buffer with current transcript lines
25 my $prevtid;
26 my ($pstart, $pend);
27 while (<>) {
28 my $curline=$_;
29 my @t=split(/\t/);
30 next unless $chr eq $t[0];
31 next if $cstrand && $cstrand ne $t[6];
32 my $tid;
33 my $gffparent;
34 if ($t[8]=~m/transcript_id\s+([^;]+)/) {
35 $tid=$1;
36 }
37 else {
38 if ($t[8]=~m/Parent=([^;]+)/) { # GFF3 format, exon/element
39 # assumes these are given in order
40 $tid=$1;
41 }
42 elsif ($t[8]=~m/ID=([^;]+)/) { # GFF3, new parent
43 $tid=$1;
44 $gffparent=1;
45 #($pstart, $pend)=($t[3],$t[4]);
46 }
47 }
48 if ($prevtid && $prevtid ne $tid) {
49 #transcript change
50 unless ($cstart>$pend || $pstart>$cend) {
51 # range overlap
52 print join('',@gffbuf);
53 }
54 @gffbuf=();
55 ($pstart, $pend) = $gffparent ? ($t[3],$t[4]) : (0,0);
56 }
57 $prevtid=$tid;
58 push(@gffbuf, $curline);
59 ($t[3], $t[4])=($t[4], $t[3]) if $t[3]>$t[4];
60 $pstart=$t[3] if ($t[3]<$pstart || $pstart==0);
61 $pend=$t[4] if $t[4]>$pend;
62 }
63 # check for last transcript read
64 if ($pstart && @gffbuf>0 && $cstart<=$pend && $pstart<=$cend) {
65 print join('',@gffbuf);
66 }

Properties

Name Value
svn:executable *