ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/gclib/scripts/gff2bed
(Generate patch)
# Line 15 | Line 15
15    select(OUTF);
16    }
17   my $tcolor= $Getopt::Std::opt_c || 0;
18 + my $tmaxscore = 0;
19   #$tcolor=hex($tcolor) if ($tcolor=~m/^0x/ || $tcolor=~m/[a-f,A-F]/);
20   # ---------------------  0     1         2         3       4      5       6       7        8  
21   my %recs; # recID =>  [ chr, strand, feat_type,  gname, tdescr, fstart, fend, [@exons], [@cds] ]
# Line 28 | Line 29
29     next if $f eq 'gene' || $f eq 'locus'; # Warning: skipping any 'gene' or 'locus' features, unconditionally
30     my $gff3_ID;
31     my $gff3_Parent;
32 +   my $score;
33     ($fstart, $fend)=($fend, $fstart) if $fend<$fstart;
34 +   ($score)=($lnum=~m/\bcov[= ]\"?([\d\.]+)/i); #or we can use fpkm of transcripts
35 +   if ($score>$tmaxscore) {
36 +      $tmaxscore=$score;
37 +      }
38     ($gff3_ID)=($lnum=~m/\bID=([^;]+)/);
39     ($gff3_Parent)=($lnum=~m/\bParent=([^;]+)/);
40     if ($gff3_ID || $gff3_Parent) { # GFF format
# Line 56 | Line 62
62             }
63           die("Error: duplicate feature $gff3_ID on $chr\n") if (exists($recs{"$chr|$gff3_ID"}));
64           my $recID="$chr|$gff3_ID";
65 <         $recs{$recID} = [$chr, $strand, $f, $gname, $tdescr, $fstart, $fend, [], [] ];
65 >         $recs{$recID} = [$chr, $strand, $f, $gname, $tdescr, $fstart, $fend, [], [], $score ];
66           next;
67           } # parent/top-level feature
68        } #GFF
# Line 95 | Line 101
101       $$ld[6]=$fend if $fend>$lend;
102       push(@{$$ld[$i]}, [$fstart, $fend, $fscore]);
103       }
104 <    else { # first time seeing this locus/gene
104 >    else { # first time seeing this transcript
105       $recs{$recID} = ($f eq 'CDS') ?
106 <           [$chr, $strand, $f, $gname, $tdescr, $fstart, $fend, [], [[$fstart, $fend, $fscore]] ] :
107 <           [$chr, $strand, $f, $gname, $tdescr, $fstart, $fend, [[$fstart, $fend, $fscore]], [] ] ;
106 >           [$chr, $strand, $f, $gname, $tdescr, $fstart, $fend, [], [[$fstart, $fend, $fscore]], $score ] :
107 >           [$chr, $strand, $f, $gname, $tdescr, $fstart, $fend, [[$fstart, $fend, $fscore]], [], $score ] ;
108              # 0       1     2     3       4        5       6     7 (exons)                8 (CDS)
109       ($gname,$tdescr)=();
110       }
# Line 135 | Line 141
141      @recs_keys=keys(%recs);
142      $rlist=\@recs_keys;
143      }
144 <    
144 >
145   foreach my $gffid (@$rlist) {
146      my $ld=$recs{$gffid};
147     # my $ld=$recs{$l} || die ("Error: locus $l found in list but not in hash!\n");
148 <   #     0       1       2       3       4      5         6    7    8
149 <   my ($chr, $strand, $ftype, $gname, $descr, $lstart, $lend, $er, $cr) = @$ld;
148 >   #     0       1       2       3       4      5         6    7    8     9
149 >   my ($chr, $strand, $ftype, $gname, $descr, $lstart, $lend, $er, $cr, $score) = @$ld;
150     my ($mstart,$mend)=($lstart, $lend);
151     my $CDexons=0;
152     my @ex;
# Line 163 | Line 169
169       push(@bstarts, $$ed[0]-$mstart);
170       push(@blens, $$ed[1]-$$ed[0]+1);
171       }
172 <   print join("\t", $chr, $mstart-1, $mend, "$gffid|$gname", 666, $strand, $cds_start-1, $cds_end, $tcolor,
172 >   if ($score) {
173 >      $tcolor=score2color($score);
174 >      
175 >      }
176 >    else {
177 >      $score=666;
178 >      }
179 >   print join("\t", $chr, $mstart-1, $mend, "$gffid|$gname", $score, $strand, $cds_start-1, $cds_end,
180 >       join(',',$tcolor,$tcolor,$tcolor),
181         scalar(@bstarts), join(',', @blens), join(',', @bstarts))."\n";
182    } #for each stored transcript
183   }
184 +
185 + sub score2color {
186 + my $c=int(255-(255.0 * ($_[0]/$tmaxscore)));
187 + $c=0 if $c<0;
188 + $c=128 if $c>=128;
189 + return $c;
190 + }

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines