Permalink
Browse files

fix merge conflict

  • Loading branch information...
2 parents 3fc8374 + 6592b11 commit a309d796dac172e0b13c7f3726d9f6176ed6f856 @lstein lstein committed Jan 27, 2013
View
6 .gitignore
@@ -0,0 +1,6 @@
+*~
+_build
+blib
+MYMETA.yml
+Build
+MYMETA.json
View
21 Changes
@@ -1,8 +1,29 @@
Revision history for Perl extension Bio::Graphics.
+2.31 Tue Sep 25 22:39:43 EDT 2012
+ - Fix infinite loop when drawing wiggle_xyplots with z-score scaling
+ over regions that do not vary much.
+
+2.30
+ - Added the glyph for FlyBase's "stacked wiggle" plot (fb_shmiggle.pm) as well as a utility script for preparing the data (index_cov_files.pl).
+
+2.29
+ - Fixes for Bio::Graphics::Wiggle to prevent crash when the "bottom" keystyle selected.
+ - Add ability to obtain Bio::Graphics::Wiggle features directly from the
+ Bio::Graphics::Wiggle::Loader without creating an intermediate gff3 file.
+
+2.28
+ - Fixes to restore density plot functionality in hybrid_plot and vista_plot.
+ - Fix handling of Bio::Graphics::Wiggle files in xyplots.
+
+2.27
+ - Fix handling of code ref subs to allow for option=>\&Packagename::functionname.
+
2.26
- Support for normalizing quantitative plots across entire track.
- Support for transparency within tracks, allowing features to overlap.
+ - Support for -color_series and -color_cycle options (see Bio::Graphics::Panel)
+ which cycle dynamically among a fixed series of background colors.
2.25
- Deprecate xyplot "histogram" subtype and default to "boxes".
View
14 MANIFEST 100755 → 100644
@@ -34,6 +34,7 @@ lib/Bio/Graphics/Glyph/ellipse.pm
lib/Bio/Graphics/Glyph/ex.pm
lib/Bio/Graphics/Glyph/extending_arrow.pm
lib/Bio/Graphics/Glyph/Factory.pm
+lib/Bio/Graphics/Glyph/fb_shmiggle.pm
lib/Bio/Graphics/Glyph/fixedwidth.pm
lib/Bio/Graphics/Glyph/flag.pm
lib/Bio/Graphics/Glyph/gene.pm
@@ -59,6 +60,7 @@ lib/Bio/Graphics/Glyph/pairplot.pm
lib/Bio/Graphics/Glyph/pentagram.pm
lib/Bio/Graphics/Glyph/phylo_align.pm
lib/Bio/Graphics/Glyph/pinsertion.pm
+lib/Bio/Graphics/Glyph/point_glyph.pm
lib/Bio/Graphics/Glyph/primers.pm
lib/Bio/Graphics/Glyph/processed_transcript.pm
lib/Bio/Graphics/Glyph/protein.pm
@@ -71,6 +73,7 @@ lib/Bio/Graphics/Glyph/repeating_shape.pm
lib/Bio/Graphics/Glyph/rndrect.pm
lib/Bio/Graphics/Glyph/ruler_arrow.pm
lib/Bio/Graphics/Glyph/saw_teeth.pm
+lib/Bio/Graphics/Glyph/scale.pm
lib/Bio/Graphics/Glyph/segmented_keyglyph.pm
lib/Bio/Graphics/Glyph/segments.pm
lib/Bio/Graphics/Glyph/smoothing.pm
@@ -96,8 +99,8 @@ lib/Bio/Graphics/Glyph/wave.pm
lib/Bio/Graphics/Glyph/weighted_arrow.pm
lib/Bio/Graphics/Glyph/whiskerplot.pm
lib/Bio/Graphics/Glyph/wiggle_box.pm
-lib/Bio/Graphics/Glyph/wiggle_density.pm
lib/Bio/Graphics/Glyph/wiggle_data.pm
+lib/Bio/Graphics/Glyph/wiggle_density.pm
lib/Bio/Graphics/Glyph/wiggle_whiskers.pm
lib/Bio/Graphics/Glyph/wiggle_xyplot.pm
lib/Bio/Graphics/Glyph/xyplot.pm
@@ -111,13 +114,15 @@ lib/Bio/Graphics/Wiggle/Loader.pm
Makefile.PL
MANIFEST
MANIFEST.SKIP
+META.json
META.yml
README
README.bioperl
scripts/contig_draw.pl
scripts/feature_draw.pl
scripts/frend.pl
scripts/glyph_help.pl
+scripts/index_cov_files.pl
scripts/render_msa.pl
scripts/search_overview.pl
t/BioGraphics.t
@@ -129,6 +134,8 @@ t/data/t1/version11.gif
t/data/t1/version11.png
t/data/t1/version12.gif
t/data/t1/version12.png
+t/data/t1/version13.gif
+t/data/t1/version13.png
t/data/t1/version2.gif
t/data/t1/version2.png
t/data/t1/version3.gif
@@ -155,6 +162,8 @@ t/data/t2/version17.gif
t/data/t2/version17.png
t/data/t2/version18.gif
t/data/t2/version18.png
+t/data/t2/version19.gif
+t/data/t2/version19.png
t/data/t2/version2.gif
t/data/t2/version2.png
t/data/t2/version3.gif
@@ -175,6 +184,8 @@ t/data/t3/version12.gif
t/data/t3/version12.png
t/data/t3/version13.gif
t/data/t3/version13.png
+t/data/t3/version14.gif
+t/data/t3/version14.png
t/data/t3/version2.gif
t/data/t3/version2.png
t/data/t3/version3.gif
@@ -188,4 +199,3 @@ t/data/t3/version8.png
t/data/t3/version9.png
t/data/wig_data.wig
t/Wiggle.t
-META.json
View
12 MANIFEST.SKIP
@@ -0,0 +1,12 @@
+^MYMETA.yml$
+^MYMETA\.json$
+^_build
+\.tar\.gz$
+^blib/
+^\.
+Build$
+^[^/]+\.pl$
+^[^/]+\.t$
+^[^/]+\.gff3
+MANIFEST\.bak
+~$
View
2 lib/Bio/Graphics.pm
@@ -2,7 +2,7 @@ package Bio::Graphics;
use strict;
use Bio::Graphics::Panel;
-our $VERSION = '2.26';
+our $VERSION = '2.31';
1;
View
24 lib/Bio/Graphics/FeatureFile.pm
@@ -434,13 +434,16 @@ END
(my $name = $args{-file}) =~ s!/!_!g;
my $cachefile = $self->cachefile($name);
if (-e $cachefile && (stat(_))[9] >= $self->file_mtime($args{-file})) { # cache is valid
+# if (-e $cachefile && -M $cachefile < 0) { # cache is valid
my $parsed_file = lock_retrieve($cachefile);
$parsed_file->initialize_code if $parsed_file->safe;
return $parsed_file;
} else {
mkpath(dirname($cachefile));
my $parsed = $self->_new(@_);
+ $parsed->initialize_code();
eval {lock_store($parsed,$cachefile)};
+ warn $@ if $@;
return $parsed;
}
@@ -1166,16 +1169,17 @@ sub code_setting {
my $setting = $self->_setting($section=>$option);
return unless defined $setting;
return $setting if ref($setting) eq 'CODE';
- if ($setting =~ /^\\&([\w:]+::)*(\w+)/) { # coderef in string form
- my $package = $1;
- my $subroutine_name = $2;
- $package ||= ($self->base2package.'::');
- my $codestring = "\\&${package}${subroutine_name}";
- my $coderef = eval $codestring;
- $self->_callback_complain($section,$option) if $@;
- $self->set($section,$option,$coderef);
- $self->set_callback_source($section,$option,$setting);
- return $coderef;
+ if ($setting =~ /^\\&([:\w]+)/) { # coderef in string form
+ my $subroutine_name = $1;
+ my $package = $self->base2package;
+ my $codestring = $subroutine_name =~ /::/
+ ? "\\&$subroutine_name"
+ : "\\&${package}\:\:${subroutine_name}" ;
+ my $coderef = eval $codestring;
+ $self->_callback_complain($section,$option) if $@;
+ $self->set($section,$option,$coderef);
+ $self->set_callback_source($section,$option,$setting);
+ return $coderef;
}
elsif ($setting =~ /^sub\s*(\(\$\$\))*\s*\{/) {
my $package = $self->base2package;
View
2 lib/Bio/Graphics/Glyph.pm
@@ -849,7 +849,7 @@ sub layout_sort {
sub layout {
my $self = shift;
return $self->{layout_height} if exists $self->{layout_height};
-
+
my @parts = $self->parts;
return $self->{layout_height} =
$self->height + $self->pad_top + $self->pad_bottom unless @parts;
View
3 lib/Bio/Graphics/Glyph/diamond.pm
@@ -2,7 +2,7 @@ package Bio::Graphics::Glyph::diamond;
# DAS-compatible package to use for drawing a colored diamond
use strict;
-use base qw(Bio::Graphics::Glyph::generic);
+use base qw(Bio::Graphics::Glyph::point_glyph);
sub my_description {
return <<END;
@@ -70,6 +70,7 @@ sub draw_component {
}
}
+
1;
__END__
View
2 lib/Bio/Graphics/Glyph/dot.pm
@@ -2,7 +2,7 @@ package Bio::Graphics::Glyph::dot;
# DAS-compatible package to use for drawing a ring or filled circle
use strict;
-use base qw(Bio::Graphics::Glyph::generic);
+use base qw(Bio::Graphics::Glyph::point_glyph);
use constant PI => 3.14159;
sub my_description {
return <<END;
View
645 lib/Bio/Graphics/Glyph/fb_shmiggle.pm
@@ -0,0 +1,645 @@
+package Bio::Graphics::Glyph::fb_shmiggle;
+
+# [ to see proper formatting set tab==2 ]
+#
+# 2009-2010 Victor Strelets, FlyBase.org
+
+use constant DEBUG => 0;
+
+#use strict;
+use GD;
+use vars '@ISA';
+use Bio::Graphics::Glyph::generic;
+BEGIN{ @ISA = 'Bio::Graphics::Glyph::generic';
+ local($colors_selected,$black,$red,$white,$grey);
+ @colors= ();
+ @colcycle= (
+ # red/orange
+ 'A62A2A', # +
+ 'CC3299', # +
+ 'FF7F00', # +
+ # yellow/brown
+ 'B87333', # +
+ # greenish
+ '6B8E23', # +
+ '238E68', # +
+ # blue/violet
+ '6982CA', # +
+ '2222AD', # +
+ );
+ %Indices= ();
+ local(*DATF);
+ }
+
+#--------------------------
+sub draw {
+ my $self = shift;
+ my $gd = shift;
+ my ($left,$top,$partno,$total_parts) = @_;
+ my $ft= $self->feature;
+ my $ftid= $ft->{id};
+ my $pn= $self->panel;
+ my $fgc = $self->fgcolor;
+ my $bgc = $self->bgcolor;
+ my($r,$g,$b);
+
+ unless( $colors_selected ) {
+ $black= $gd->colorClosest(0,0,0);
+ $yellow= $gd->colorClosest(255,255,0);
+ $white= $gd->colorClosest(255,255,255);
+ $red= $gd->colorClosest(255,0,0);
+ $lightgrey= $gd->colorClosest(225,225,225);
+ $grey= $gd->colorClosest(200,200,200);
+ $darkgrey= $gd->colorClosest(125,125,125);
+ foreach my $ccc ( @colcycle ) {
+ if( $ccc=~/^[#]?(\S\S)(\S\S)(\S\S)$/ ) {
+ ($r,$g,$b)= ($1,$2,$3);
+ eval('$r=hex("0x'.$r.'");');
+ eval('$g=hex("0x'.$g.'");');
+ eval('$b=hex("0x'.$b.'");');
+ }
+ my $c= $gd->colorClosest($r,$g,$b);
+ push(@colors,$c);
+ }
+ $colors_selected= 1;
+ }
+
+ my($pnstart,$pnstop)= ($pn->start,$pn->end); # in seq coordinates
+ my $nseqpoints= $pnstop - $pnstart + 1;
+ my($xf1,$yf1,$xf2,$yf2)= $self->calculate_boundaries($left,$top);
+ my $leftpad= $pn->pad_left;
+ my $datadir= $ENV{SERVER_PATH} . $ft->{datadir};
+
+ my($start,$stop)= $pn->location2pixel(($ft->{start},$ft->{end}));
+ my $ftscrstop= $stop + $leftpad;
+ my $ftscrstart= $start + $leftpad;
+
+ my $chromosome= $ft->{ref};
+ #warn("pn start stop leftpad nseq dir = $pnstart $pnstop $leftpad $datadir $chromosome\n");
+ my $flipped= $self->{flip} ? 1 : 0;
+ my($subsets,$subsetsnames,$signals)= $self->getData($ft,$datadir,$chromosome,$pnstart,$pnstop,$xf1,$xf2,$flipped);
+
+ my $poly_pkg = $self->polygon_package;
+
+ my @orderedsubsets= @{$subsets};
+ my $nsets= $#orderedsubsets+1;
+ my($xstep,$ystep)= (2,int(100.0/$nsets));
+ $ystep= 7 unless $ystep>=7; # empiricaly found - to read lines of tiny fonts
+ $ystep= 12 if $ystep>12; # empirically found - to preserve topo feel when number of subsets is small
+ $ystep= 7 if $ystep>7; # tmp unification
+ my($xw,$yw)= ( $nsets*$xstep, ($nsets-1)*$ystep );
+
+ my $polybg= $poly_pkg->new();
+ $polybg->addPt($xf1,$yf2-$yw);
+ $polybg->addPt($xf2,$yf2-$yw);
+ $polybg->addPt($xf2-$xw, $yf2);
+ $polybg->addPt($xf1-$xw, $yf2);
+ $gd->filledPolygon($polybg,$lightgrey); # background
+ for( my $xx= $xf1+2; $xx<$xf2; $xx+=6 ) { $gd->line($xx,$yf2-$yw,$xx-$xw,$yf2,$grey); } # grid-helper
+
+ my $xshift= 0;
+ my $yshift= $nsets * $ystep;
+ ($r,$g,$b)= (10,150,80);
+ my $colcycler= 0;
+ my @screencoords= @{$signals->{screencoords}};
+ my $max_signal= 30;
+ my $koeff= 4;
+ if( exists $signals->{max_signal} ) {
+ $max_signal= $signals->{max_signal};
+ $koeff= 80.0/$max_signal;
+ }
+ my $predictor_cutoff= int($max_signal*0.95); # empirically found
+ my @prevx= ();
+ my @prevy= ();
+ my @prevvals= ();
+ my $profilen= 0;
+ my %SPEEDUP= ();
+ foreach my $subset ( @orderedsubsets ) {
+ my $edgecolor= ($xshift==0) ? $black : $yellow;
+ #my $color= ($xshift==0) ? $darkgrey : $gd->colorClosest(oct($r),oct($g),oct($b));
+ my $color= ($profilen==0) ? $darkgrey : $colors[$colcycler];
+ $xshift -= $xstep;
+ $yshift -= $ystep;
+ my @values= @{$signals->{$subset}};
+ my($xold,$yold)= ($xf1+$xshift,$yf2-$yshift+1);
+ my $xpos= 0;
+ my $poly= $poly_pkg->new();
+ $poly->addPt($xold,$yold+1);
+ my @allx= ($xold);
+ my @ally= ($yold);
+ my @allvals= (0);
+ my $runx= $xf1 + $xshift;
+ foreach my $val ( @values ) {
+ $scrx += $leftpad;
+ my $x= $screencoords[$xpos] + $xshift;
+ my $visval;
+ if( exists $SPEEDUP{$val} ) { $visval= $SPEEDUP{$val}; }
+ else { $visval= int($val*$koeff); $SPEEDUP{$val}= $visval; }
+ my $y= $yf2 - $yshift - $visval;
+ push(@allx,$x);
+ push(@ally,$y);
+ push(@allvals,$visval);
+ if( $xpos>0 ) {
+ $poly->addPt($x,$y+1);
+ }
+ ($xold,$yold)= ($x,$y);
+ $xpos++;
+ }
+ $poly->addPt($xf2+$xshift, $yf2-$yshift+1);
+ $gd->filledPolygon($poly,$color) unless $profilen==0; # not on MAX predictor
+ ($xold,$yold)= ($allx[0],$ally[0]);
+ for( my $en=1; $en<=$#allx; $en++ ) {
+ my $x= $allx[$en];
+ my $y= $ally[$en];
+ $gd->line($xold,$yold,$x,$y,$edgecolor);
+ ($xold,$yold)= ($x,$y);
+ }
+ if( $profilen==0 ) { # drawing mRNA (cutoff-based) predictor on MAX subset
+ my($xxx,$yyy)= ($allx[1]-1,$yf2-$yw);
+ $gd->line($xxx-4,$yyy,$xxx-2,$yyy,$black);
+ $gd->string(GD::Font->Tiny,$xxx-12, $yyy-3,'0',$black);
+ $gd->line($xxx-2,$yyy,$xxx-2,$yyy-50,$black);
+ $gd->line($xxx-4,$yyy-47,$xxx-2,$yyy-50,$black);
+ $gd->line($xxx,$yyy-47,$xxx-2,$yyy-50,$black);
+ $gd->line($xxx-4,$yyy-44,$xxx-2,$yyy-44,$black);
+ $gd->string(GD::Font->Tiny,$xxx-18, $yyy-47,$max_signal,$black);
+ my($inexon,$exstart,$exend,$ymax)= (0,0,0,999);
+ for( my $en=1; $en<=$#allx; $en++ ) {
+ my $y= $ally[$en];
+ $ymax= $y if $y < $ymax;
+ if( $allvals[$en]>=$predictor_cutoff ) {
+ my $x= $allx[$en];
+ unless( $inexon ) { $inexon= 1; $exstart= $x; } # start exon
+ $exend= $x;
+ }
+ elsif( $inexon ) { # end exon and draw it
+ $inexon= 0;
+ $ymax -= 6;
+ my $allowedymax= $yf2-$yshift-45;
+ $ymax= $allowedymax if $ymax < $allowedymax; # set limit for huge peaks
+ $gd->line($exstart,$ymax,$exstart,$ymax+2,$red);
+ $gd->line($exstart,$ymax,$exend,$ymax,$red);
+ $gd->line($exend,$ymax,$exend,$ymax+2,$red);
+ ($inexon,$exstart,$exend,$ymax)= (0,0,0,999);
+ }
+ }
+ if( $inexon ) { # exon which ends beyond this screen
+ $ymax -= 6;
+ my $allowedymax= $yf2-$yshift-45;
+ $ymax= $allowedymax if $ymax < $allowedymax; # set limit for huge peaks
+ $gd->line($exstart,$ymax,$exstart,$ymax+2,$red);
+ $gd->line($exstart,$ymax,$exend,$ymax,$red);
+ }
+ }
+ if( 0 && $profilen>1 ) { # blocked - drawing roof lines doesn't work well with this view..
+ for( my $en=3; $en<=$#allx; $en+=6 ) {
+ next if $allvals[$en]==0 || $prevvals[$en]==0;
+ my $y= $ally[$en];
+ $yold= $prevy[$en];
+ if( $yold<$y ) {
+ my $x= $allx[$en];
+ $xold= $prevx[$en];
+ $gd->line($xold,$yold,$x,$y,$darkgrey);
+ }
+ }
+ }
+ $gd->string(GD::Font->Tiny,$xf2+$xshift+3, $yf2-$yshift-5,$subsetsnames->{$subset},$color);
+ $colcycler++;
+ $colcycler= 0 if $colcycler>$#colors;
+ unless( $profilen==0 ) { @prevx= @allx; @prevy= @ally; @prevvals= @allvals; }
+ $profilen++;
+ }
+
+ return;
+}
+
+#--------------------------
+sub getData {
+ my $self = shift;
+ my($ft,$datadir,$chromosome,$start,$stop,$scrstart,$scrstop,$flipped) = @_;
+ my %Signals= ();
+ $self->openDataFiles($datadir);
+ my @subsets= (exists $ft->{'subsetsorder'}) ? @{$ft->{'subsetsorder'}} : sort split(/\t+/,$Indices{'subsets'});
+ shift(@subsets) if $subsets[0] eq 'MAX';
+ warn("subsets: @subsets\n") if DEBUG;
+ my %SubsetsNames= (exists $ft->{'subsetsnames'}) ? %{$ft->{'subsetsnames'}} : map { $_, $_ } @subsets;
+ $SubsetsNames{MAX}= 'MAX';
+ my $screenstep= ($scrstop-$scrstart+1) * 1.0 / ($stop-$start+1);
+ my $donecoords= 0;
+ foreach my $subset ( @subsets ) {
+ my $nstrings= 0;
+ # scan seq ranges offsets to see where to start reading
+ my $key= $subset.':'.$chromosome;
+ my $poskey= $key.':offsets';
+ my $ranges_pos= (exists $Indices{$poskey}) ? int($Indices{$poskey}) : -1;
+ if( $ranges_pos == -1 ) { next; } # no such signal..
+ warn(" positioning for $poskey starts at $ranges_pos\n") if DEBUG;
+ if( $start>=1000000 ) {
+ my $bigstep= int($start/1000000.0);
+ if( exists $Indices{$key.':offsets:'.$bigstep} ) {
+ my $jumpval= $Indices{$key.':offsets:'.$bigstep};
+ warn(" jump in offset search to $jumpval\n") if DEBUG;
+ $ranges_pos= int($jumpval); }
+ }
+ seek(DATF,$ranges_pos,0);
+ my($offset,$offset1)= (0,0);
+ my $lastseqloc= -999999999;
+ my $useoffset= 0;
+ while( (my $strs=<DATF>) ) {
+ $nstrings++ if DEBUG;
+ if( DEBUG ) {
+ chop($strs); warn(" positioning read for coord $start ($strs)\n"); }
+ last unless $strs=~m/^(-?\d+)[ \t]+(\d+)/;
+ my($seqloc,$fileoffset)= ($1,$2);
+ if( DEBUG ) {
+ chop($strs); warn(" positioning read for $poskey => $seqloc, $fileoffset ($strs)\n"); }
+ $offset1= $offset;
+ $offset= $fileoffset;
+ $lastseqloc= $seqloc;
+ if( $seqloc > $start ) { $useoffset= int($offset1); last; }
+ }
+ warn(" will use offset $useoffset\n") if DEBUG;
+ warn(" (scanned $nstrings offset strings)\n") if DEBUG;
+ if( $useoffset==0 ) { # data offset cannot be 0 - means didn't find where to read required data..
+ next;
+ my @emptyvals= ();
+ for( my $ii= $scrstart; $ii++ <= $scrstop; ) { push(@emptyvals,0); }
+ $Signals{$subset}= \@emptyvals;
+ }
+ $nstrings= 0;
+ # read signal profile
+ seek(DATF,$useoffset,0);
+ $lastseqloc= -999999999;
+ my $lastsignal= 0;
+ my($scrx,$scrxold)= ($scrstart,$scrstart-1);
+ my $runmax= 0;
+ my @values= ();
+ my @xscreencoords= ();
+ while( (my $str=<DATF>) ) {
+ $nstrings++ if DEBUG;
+ unless( $str=~m/^(-?\d+)[ \t]+(\d+)/ ) {
+ warn(" header read: $str") if DEBUG;
+ last; # because no headers were indexed at the beginning of data packs
+ }
+ my($seqloc,$signal)= ($1,$2);
+ warn(" signal read: $seqloc, $signal line: $str") if DEBUG;
+ last if $lastseqloc > $seqloc; # just in case, as all sits merged in one file..
+ if( $seqloc>=$start ) { # current is the next one after the one we need to start from..
+ unless( $lastseqloc== -999999999 ) { # expand previous
+ $lastseqloc= $start-2 if $lastseqloc<$start; # limit empty steps (they may start from -200000)
+ while( $lastseqloc < $seqloc ) { # until another (one we just retrieved) wiggle reading
+ last if $lastseqloc > $stop; # end of subset data
+ next if $lastseqloc++ < $start;
+ # we have actual new seq position in our required range
+ my $scrpos= int($scrx);
+ $runmax= $lastsignal if $runmax < $lastsignal;
+ if( $scrpos != $scrxold ) { # we have actual new seq _and_ screen position
+ push(@values,$runmax);
+ push(@xscreencoords,$scrpos) unless $donecoords;
+ $scrxold= $scrpos;
+ $runmax= 0;
+ }
+ $scrx += $screenstep; # remember - it is not integer
+ }
+ }
+ }
+ ($lastseqloc,$lastsignal)= ($seqloc,$signal);
+ last if $seqloc > $stop; # end of subset data
+ }
+ if( $lastseqloc < $stop ) { # if on the end of signal profile, but still in screen range
+ # just assume that we are getting one more reading with signal == 0
+ my $signal= 0;
+ while( $lastseqloc++ < $stop ) {
+ my $scrpos= int($scrx);
+ if( $scrpos != $scrxold ) { # we have actual new seq _and_ screen position
+ push(@values,$signal);
+ push(@xscreencoords,$scrpos) unless $donecoords;
+ $scrxold= $scrpos;
+ }
+ $scrx += $screenstep;
+ }
+ }
+ warn(" (scanned $nstrings signal strings)\n") if DEBUG;
+ $nstrings= 0;
+ if( $flipped ) {
+ my @ch= reverse @values; @values= @ch;
+ }
+ warn(" ".$subset."=> ".@values." values @values\n") if DEBUG && $#values<1000;
+ $Signals{$subset}= \@values;
+ $Signals{screencoords}= \@xscreencoords unless $donecoords;
+ $donecoords= 1;
+ } # foreach my $subset ( @subsets ) {
+ if( exists $Indices{max_signal} ) {
+ $Signals{max_signal}= $Indices{max_signal};
+ warn(" max_signal=> ".$Indices{max_signal}." \n") if DEBUG;
+ }
+ # prepare MAX profile - will be used as a base for exon/UTR prediction
+ my @maxprofile= ();
+ my @ruler= @{$Signals{screencoords}};
+ for( my $npos= 0; $npos<=$#ruler; $npos++ ) {
+ my $maxval= 0;
+ foreach my $subset ( @subsets ) {
+ my $p= $Signals{$subset};
+ my $val= $p->[$npos];
+ $maxval= $val if $maxval < $val;
+ }
+ push(@maxprofile,$maxval);
+ }
+ $Signals{MAX}= \@maxprofile;
+ warn(" MAX=> ".@maxprofile." values @maxprofile\n") if DEBUG && $#maxprofile<1000;
+ unshift(@subsets,'MAX');
+
+ return(\@subsets,\%SubsetsNames, \%Signals);
+}
+
+#--------------------------
+sub openDataFiles {
+ my $self = shift;
+ my $datadir= shift;
+ $datadir.= '/' unless $datadir=~m|/$|;
+ my $datafile= $datadir.'data.cat';
+ open(DATF,$datafile) || warn("cannot open $datafile\n");
+ use BerkeleyDB; # caller should already used proper 'use lib' command with path
+ my $bdbfile= $datadir . 'index.bdbhash';
+ tie %Indices, "BerkeleyDB::Hash", -Filename => $bdbfile, -Flags => DB_RDONLY || warn("can't read BDBHash $bdbfile\n");
+ if( DEBUG ) { foreach my $kk ( sort keys %Indices ) { warn(" $kk => ".$Indices{$kk}."\n"); } }
+ return;
+}
+
+1;
+
+
+=pod
+
+=head1 TopoView Glyph
+
+=begin html
+
+<i>Warning:</i> This software is still in the developmental stage and is distributed
+"as is", without packaging and in the same exact condition as currently used by FlyBase.
+You are free to use it, modify and develop further, but proper reference
+to the original author and FlyBase is required.
+
+<p>
+
+"fb_shmiggle.pm" TopoView (AKA shmiggle) glyph was developed for fast
+3D-like demonstration of RNA-seq data consisting of multiple
+individual subsets. Main purposes were to compact presentation as
+much as possible (in one reasonably sized track) and
+to allow easy visual detection of coordinated behavior
+of the expression profiles of different subsets.
+
+<p>
+
+It was found that log2 conversion dramatically changes
+perception of expression profiles and kind of illuminates
+coordinated behavior of different subsets. Glyph and data
+indexer/formatter were in fact modified with the assumption
+that final data produced by indexer/formatter will always
+be a log2 conversion of the original coverage, therefore
+represented by short integer with values in range of 0-200
+or so.
+
+<p>
+
+Comparing performance (retrieval of several Kbp of data profiles
+for several subsets of some RNA-seq experiment) of wiggle binary
+method and of several possible alternatives, it was discovered that
+one of the approaches remarkably outperforms wiggle bin method
+(although it requires several times more space for formatted data
+storage). Optimal storage/retrieval method stores all experiment
+data (all subsets of the experiment) in one text file, where
+structure of the file in fact is one of the most simple wiggle
+(coverage files) formats with the addition of some positioning
+data (two-column format, without runlength specification, without
+omission of zero values). This is the only format which glyph is able
+to handle (there are many reasons for that) so any modification
+of indexer/formatter _must_ produce exact equivalent of that
+format. In my experience, 90% of the debugging with new incoming
+data was related to
+the problems of that exact format conversion. Example of the formatted
+data:
+
+<br>
+
+<pre>
+# subset=BS107_all_unique chromosome=2LHet
+-200000 0
+0 0
+19955 1
+19959 0
+19967 2
+19972 0
+19977 2
+20027 0
+20031 2
+20035 0
+20043 1
+20045 0
+20049 1
+20055 0
+20062 2
+20069 0
+20073 2
+20082 0
+20097 3
+20115 0
+20125 3
+20127 0
+20134 3
+20139 0
+20140 3
+20144 0
+20145 3
+20150 0
+20157 3
+20162 0
+20172 3
+20183 0
+</pre>
+
+<p>
+
+Glyph is supplied with a "index_cov_files.pl" data indexer/formatter
+which is converting original coverage (wiggle) files into data structure which will
+be used for fast retrieval. You should run this script in some separate directory,
+containing original coverage files (gzipped form works too). After it finishes,
+directory will contain two new files: data.cat and index.bdbhash. Both files required
+for data retrieval by glyph. Files can be moved freely between different directories
+or even operational systems (Mac and PC included, I think). Content of the dat file
+is subject of accurate check - this is if you want to avoid long debugging sessions
+on the level of running GBrowse. Size of files is quite big, but in my experience it
+is like twice less than gzipped size of all initial coverage files - which is quite
+acceptable.
+
+<p>
+
+Example of GBrowse conf file insert (shows actual FlyBase config sections for
+Baylor and modENCODE RNA-seq tracks):
+
+<br>
+<pre>
+[baylor_wiggle]
+feature = RNAseq_profile:Baylor
+glyph = fb_shmiggle
+height = 124
+bgcolor = sub { my $f= shift;
+ $f->{datadir}= '/.data/genomes/dmel/current/rnaseq-gff/baylor/'; # trick it this way..
+ my @subsetsorder= qw(
+ E2-4hr
+ E2-16hr
+ E2-16hr100
+ E14-16hr
+ L
+ L3i
+ L3i100
+ P
+ P3d
+ MA3d
+ FA3d
+ A17d
+ );
+ $f->{subsetsorder}= \@subsetsorder;
+ return 'lightgrey';
+ }
+key = Baylor group RNA-seq coverage by subsets (devel.stages) [log2 converted]
+category = RNA-seq data
+label = ""
+title = ""
+link = sub { my $f= shift;
+ my $id= $f->{'id'};
+ my $lnk="javascript:void(0);";
+ "$lnk\" id=\"$id\" onmouseover=\"showdata_description('Baylor');return false;\" onmouseout=\"delsumm_overlib();";
+ }
+
+[celniker_wiggle]
+feature = RNAseq_profile:Celniker
+glyph = fb_shmiggle
+height = 250
+bgcolor = sub { my $f= shift;
+ $f->{datadir}= '/.data/genomes/dmel/current/rnaseq-gff/celniker/'; # trick it this way..
+ my @subsetsorder= qw(
+ BS40_all_unique
+ BS43_all_unique
+ BS46_all_unique
+ BS49_all_unique
+ BS54_all_unique
+ BS55_all_unique
+ BS58_all_unique
+ BS62_all_unique
+ BS66_all_unique
+ BS67_all_unique
+ BS71_all_unique
+ BS73_all_unique
+ BS107_all_unique
+ BS111_all_unique
+ BS113_all_unique
+ BS196_all_unique
+ BS200_all_unique
+ BS203_all_unique
+ BS129_all_unique
+ BS133_all_unique
+ BS136_all_unique
+ BS137_all_unique
+ BS140_all_unique
+ BS143_all_unique
+ BS150_all_unique
+ BS156_all_unique
+ BS162_all_unique
+ BS153_all_unique
+ BS159_all_unique
+ BS165_all_unique
+ );
+ $f->{subsetsorder}= \@subsetsorder;
+ my %subsetsnames= qw(
+ BS40_all_unique em0-2hr
+ BS43_all_unique em2-4hr
+ BS46_all_unique em4-6hr
+ BS49_all_unique em6-8hr
+ BS54_all_unique em8-10hr
+ BS55_all_unique em10-12hr
+ BS58_all_unique em12-14hr
+ BS62_all_unique em14-16hr
+ BS66_all_unique em16-18hr
+ BS67_all_unique em18-20hr
+ BS71_all_unique em20-22hr
+ BS73_all_unique em22-24hr
+ BS107_all_unique L1
+ BS111_all_unique L2
+ BS113_all_unique L3_12hr
+ BS196_all_unique L3_PS1-2
+ BS200_all_unique L3_PS3-6
+ BS203_all_unique L3_PS7-9
+ BS129_all_unique WPP
+ BS133_all_unique WPP_12hr
+ BS136_all_unique WPP_24hr
+ BS137_all_unique WPP_2days
+ BS140_all_unique WPP_3days
+ BS143_all_unique WPP_4days
+ BS150_all_unique AdM_Ecl_1days
+ BS156_all_unique AdM_Ecl_5days
+ BS162_all_unique AdM_Ecl_30days
+ BS153_all_unique AdF_Ecl_1days
+ BS159_all_unique AdF_Ecl_5days
+ BS165_all_unique AdF_Ecl_30days
+ );
+ $f->{subsetsnames}= \%subsetsnames;
+ return 'lightgrey';
+ }
+key = modENCODE Transcription Group RNA-seq coverage (unique reads only) by subsets (devel. stages) [log2 converted]
+category = RNA-seq data
+label = ""
+title = ""
+link = sub { my $f= shift;
+ my $id= $f->{'id'};
+ my $lnk="javascript:void(0);";
+ "$lnk\" id=\"$id\" onmouseover=\"showdata_description('Celniker');return false;\" onmouseout=\"delsumm_overlib();";
+ }
+</pre>
+<br>
+
+In configuration, it is very important to set 'datadir' variable (relative
+to server DOCUMENT_ROOT) so that glyph will know where to take data and index.
+
+<br>
+<br>
+
+Setting 'subsetsorder' allows you to display expression profiles of subsets in
+some predefined order. If setting omitted, glyph will display sets in alphabetical
+order of the initial subsets names.
+
+<br>
+<br>
+
+Setting 'subsetsnames' allows to rename subsets (very important as in most cases
+workflow names of subsets are unsutable for intelligent data display to end users).
+If setting omitted, initial subsets names will be used for display.
+
+<p>
+
+For the glyph to be properly activated, you need to insert in all of your GFF files
+(ones for which you have RNA-seq data) virtual contig-long features which will activate
+expression data display. To cover whole range of the contig (chromosome arm), it is
+better to use coordinates presented in 'sequence-region' definition at the top of GFF file.
+Example of such feature lines for FlyBase data is shown below:
+
+<p>
+<pre>
+2LHet Baylor RNAseq_profile 1 368874 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+2L Baylor RNAseq_profile 1 23011544 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+2RHet Baylor RNAseq_profile 1 3288763 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+2R Baylor RNAseq_profile 1 21146708 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+3LHet Baylor RNAseq_profile 1 2555493 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+3L Baylor RNAseq_profile 1 24543557 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+3RHet Baylor RNAseq_profile 1 2517509 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+3R Baylor RNAseq_profile 1 27905053 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+4 Baylor RNAseq_profile 1 1351857 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+XHet Baylor RNAseq_profile 1 204113 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+X Baylor RNAseq_profile 1 22422827 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+YHet Baylor RNAseq_profile 1 347040 . + . Comment=This is a reference feature for RNAseq wiggle tracks
+</pre>
+<p>
+
+Questions about TopoView glyph should be directed to Victor Strelets (strelets@bio.indiana.edu).
+
+=end html
View
122 lib/Bio/Graphics/Glyph/hybrid_plot.pm
@@ -1,11 +1,10 @@
package Bio::Graphics::Glyph::hybrid_plot;
use strict;
-use base qw(Bio::Graphics::Glyph::wiggle_xyplot Bio::Graphics::Glyph::smoothing);
+use base qw(Bio::Graphics::Glyph::wiggle_xyplot);
use constant DEBUG=>0;
use constant NEGCOL=>"orange";
use constant POSCOL=>"blue";
-
our $VERSION = '1.0';
sub my_options {
@@ -72,19 +71,13 @@ sub draw {
my $height = $bottom - $top;
my $feature = $self->feature;
my $set_flip = $self->option('flip_sign') | 0;
-
-
+
#Draw individual features for reads (unlike wiggle features reads will have scores)
my $t_id = $feature->method;
if($t_id && $t_id eq $self->_check_uni){return Bio::Graphics::Glyph::generic::draw_component($self,@_);}
#Draw multiple graph if we don't have a score
- my @wiggles = ();
- foreach ('A'..'Z') {
- my $filename = 'wigfile'.$_;
- my ($wiggle) = $feature->get_tag_values('wigfile'.$_);
- push (@wiggles, $wiggle);
- }
+ my @wiggles = $self->get_wiggles($feature);
my ($fasta) = $feature->get_tag_values('fasta');
my($scale,$y_origin,$min_score,$max_score);
@@ -93,50 +86,99 @@ sub draw {
#Depending on what we have (wiggle or BigWig) pick the way to paint the signal graph
for(my $w = 0; $w < @wiggles; $w++){
- if ($w > 0) {
- $self->configure(-pos_color, NEGCOL);
- $self->configure(-no_grid, 1);
- } else {
- $self->configure(-pos_color, POSCOL);
- }
- if ($wiggles[$w] =~ /\.wi\w{1,3}$/) {
- $self->draw_wigfile($feature,$wiggles[$w],@_);
- } elsif ($wiggles[$w] =~ /\.bw$/ && $fasta) {
- my $flip = ($w > 0 && $set_flip) ? -1 : 1;
- use Bio::DB::BigWig 'binMean';
- use Bio::DB::Sam;
- my $wig = Bio::DB::BigWig->new(-bigwig => "$wiggles[$w]",
- -fasta => Bio::DB::Sam::Fai->open("$fasta"));
-
- my ($summary) = $wig->features(-seq_id => $feature->segment->ref,
- -start => $self->panel->start,
- -end => $self->panel->end,
- -type => 'summary');
-
- my $stats = $summary->statistical_summary($self->width);
- my @vals = map {$_->{validCount} ? $_->{sumData}/$_->{validCount}*$flip:0} @$stats;
- $self->draw_coverage($self,\@vals,@_);
- }
+ if ($w > 0) {
+ $self->configure('bgcolor', NEGCOL);
+ $self->configure('no_grid', 1);
+ } else {
+ $self->configure('bgcolor', POSCOL);
+ }
+ if ($wiggles[$w] =~ /\.wi\w{1,3}$/) {
+ $self->draw_wigfile($feature,$wiggles[$w],@_);
+ } elsif ($wiggles[$w] =~ /\.bw$/) {
+ my $flip = ($w > 0 && $set_flip) ? -1 : 1;
+ eval "require Bio::DB::BigWig;1" or die $@;
+ eval "require Bio::DB::Sam; 1" or die $@;
+ my @args = (-bigwig => "$wiggles[$w]");
+ push @args,(-fasta => Bio::DB::Sam::Fai->open($fasta)) if $fasta;
+ my $wig = Bio::DB::BigWig->new(@args);
+ my ($summary) = $wig->features(-seq_id => $feature->segment->ref,
+ -start => $self->panel->start,
+ -end => $self->panel->end,
+ -type => 'summary');
+ my $stats = $summary->statistical_summary($self->width);
+ my @vals = map {$_->{validCount} ? $_->{sumData}/$_->{validCount}*$flip:0} @$stats;
+ $self->_draw_coverage($summary,\@vals,@_);
+ }
}
}
-1;
+sub get_wiggles {
+ my $self = shift;
+ my $feature = shift;
+ my @wiggles;
+ foreach ('A'..'Z') {
+ my $filename = 'wigfile'.$_;
+ my ($wiggle) = $feature->get_tag_values('wigfile'.$_);
+ push (@wiggles, $wiggle) if $wiggle;
+ }
+ return @wiggles;
+}
+sub minmax {
+ my $self = shift;
+ my $parts = shift;
+
+ my $autoscale = $self->option('autoscale') || 'local';
+
+ my $min_score = $self->min_score unless $autoscale eq 'z_score';
+ my $max_score = $self->max_score unless $autoscale eq 'z_score';
+
+ my $do_min = !defined $min_score;
+ my $do_max = !defined $max_score;
+
+ my @wiggles = $self->get_wiggles($self->feature);
+ my ($min,$max,$mean,$stdev);
+ my @args = (-seq_id => (eval{$self->feature->segment->ref}||''),
+ -start => $self->panel->start,
+ -end => $self->panel->end,
+ -type => 'summary');
+
+ for my $w (@wiggles) {
+ my ($a,$b,$c,$d);
+ if ($w =~ /\.bw$/) {
+ eval "require Bio::DB::BigWig;1" or die $@;
+ my $wig = Bio::DB::BigWig->new(-bigwig=>$w) or next;
+ ($a,$b,$c,$d) = $self->bigwig_stats($autoscale,$wig->features(@args));
+ } elsif ($w =~ /\.wi\w{1,3}$/) {
+ eval "require Bio::Graphics::Wiggle;1" or die $@;
+ my $wig = Bio::Graphics::Wiggle->new($w);
+ ($a,$b,$c,$d) = $self->wig_stats($autoscale,$wig);
+ }
+ $min = $a if !defined $min || $min > $a;
+ $max = $b if !defined $max || $max < $b;
+ $mean += $c;
+ $stdev += $d**2;
+ }
+ $stdev = sqrt($stdev);
+
+ $min_score = $min if $do_min;
+ $max_score = $max if $do_max;
+ return $self->sanity_check($min_score,$max_score,$mean,$stdev);
+}
+
+1;
__END__
=head1 NAME
-
-Bio::Graphics::Glyph::hybrid_plot - An xyplot plot drawing dual graph using data from two wiggle files per track
+Bio::Graphics::Glyph::hybrid_plot - An xyplot plot drawing dual graph using data from two or more wiggle files per track
=head1 SYNOPSIS
-
See <Bio::Graphics::Panel> <Bio::Graphics::Glyph> and <Bio::Graphics::Glyph::wiggle_xyplot>.
=head1 DESCRIPTION
-
Note that for full functionality this glyph requires Bio::Graphics::Glyph::generic (generic glyph is used for drawing individual
matches for small RNA alignments at a high zoom level, specified by semantic zooming in GBrowse conf file)
Unlike the regular xyplot, this glyph draws two overlapping graphs
@@ -165,6 +207,8 @@ sequences (signal quality) aligned with the current region so the user can see
the difference between accumulated signal from overlapping multiple matches
(which may likely be just a noise from products of degradation) and high-quality
signal from unique sequences.
+
+For a third wiggle file use the attribute "wigfileC" and so forth.
It is essential that wigfile entries in gff file do not have score, because
score used to differentiate between data for dual graph and data for matches
View
17 lib/Bio/Graphics/Glyph/point_glyph.pm
@@ -0,0 +1,17 @@
+package Bio::Graphics::Glyph::point_glyph;
+
+use strict;
+use base qw(Bio::Graphics::Glyph::generic);
+
+sub box {
+ my $self = shift;
+ my @result = $self->SUPER::box();
+ return @result unless $self->option('point') && $result[2]-$result[0] < 3;
+ my $h = $self->option('height')/2;
+ my $mid = int(($result[2]+$result[0])/2);
+ $result[0] = $mid-$h-1; # fudge a little to make it easier to click on
+ $result[2] = $mid+$h+1;
+ return @result;
+}
+
+1;
View
7 lib/Bio/Graphics/Glyph/track.pm
@@ -35,9 +35,10 @@ sub draw {
$self->normalize_track(@parts);
# dynamic assignment of colors
- if ($self->option('color_series')) {
- my @color_series =
- qw(red blue green yellow orange brown aqua black fuchsia green lime maroon navy olive purple silver teal magenta);
+ if ($self->option('color_series') || $self->option('color_cycle')) {
+ my $series = $self->option('color_cycle');
+ $series ||= 'red blue green yellow orange brown aqua black fuchsia green lime maroon navy olive purple silver teal magenta';
+ my @color_series = ref($series) eq 'ARRAY' ? @$series : split /\s+/,$series;
my $index = 0;
my %color_cache;
my $closure = sub {
View
2 lib/Bio/Graphics/Glyph/triangle.pm
@@ -2,7 +2,7 @@ package Bio::Graphics::Glyph::triangle;
# DAS-compatible package to use for drawing a triangle
use strict;
-use base qw(Bio::Graphics::Glyph::generic);
+use base qw(Bio::Graphics::Glyph::point_glyph);
sub pad_left {
my $self = shift;
View
78 lib/Bio/Graphics/Glyph/vista_plot.pm
@@ -1,12 +1,12 @@
package Bio::Graphics::Glyph::vista_plot;
use strict;
-use base qw(Bio::Graphics::Glyph::wiggle_data
- Bio::Graphics::Glyph::wiggle_xyplot
+use base qw(Bio::Graphics::Glyph::wiggle_xyplot
Bio::Graphics::Glyph::wiggle_density
Bio::Graphics::Glyph::wiggle_whiskers
Bio::Graphics::Glyph::heat_map
- Bio::Graphics::Glyph::smoothing);
+ Bio::Graphics::Glyph::smoothing);
+
our $VERSION = '1.0';
@@ -46,8 +46,8 @@ sub my_options {
'vista',
"What to show, peaks or signal, both (vista plot) or density graph."],
graph_type => [
- ['whiskers','histogram','line','points','linepoints'],
- 'histogram',
+ ['whiskers','histogram','boxes','line','points','linepoints'],
+ 'boxes',
"Type of signal graph to show."],
alpha => [
'integer',
@@ -123,8 +123,10 @@ sub bigwig_summary {
sub global_mean_and_variance {
my $self = shift;
if (my $wig = $self->wig) {
- return ($wig->mean,$wig->stdev);
- } elsif (my $sum = $self->bigwig_summary){
+ my @result = eval{($wig->mean,$wig->stdev)};
+ return @result if @result;
+ }
+ if (my $sum = $self->bigwig_summary){
use Bio::DB::BigWig qw(binMean binStdev);
my $stats = $sum->statistical_summary(1);
return eval{(binMean($stats->[0]),binStdev($stats->[0]))};
@@ -158,9 +160,9 @@ sub draw {
peak => (eval{$feature->get_tag_values('peak_type')})[0],
fasta=> (eval{$feature->get_tag_values('fasta')})[0]);
$self->panel->startGroup($gd);
- $self->configure(opacity => 0.5) if $only_show eq 'vista';
+
$self->draw_signal($only_show,\%features,@_) if $only_show =~ /signal|density|vista/;
- $self->draw_peaks(\%features,@_) if $features{peak} && $only_show =~ /peaks|vista/;
+ $self->draw_peaks(\%features,@_) if $features{peak} && $only_show =~ /peaks|vista|both/;
$self->Bio::Graphics::Glyph::xyplot::draw_label(@_) if $self->option('label');
$self->draw_description(@_) if $self->option('description');
$self->panel->endGroup($gd);
@@ -204,9 +206,9 @@ sub draw_signal {
my @vals = map {$_->{validCount} ? Bio::DB::BigWig::binMean($_) : 0} @$stats;
$self->bigwig_summary($summary);
if ($signal_type eq 'density') {
- $self->Bio::Graphics::Glyph::wiggle_density::draw_coverage($summary,\@vals,@_);
+ $self->Bio::Graphics::Glyph::wiggle_density::_draw_coverage($summary,\@vals,@_);
} else {
- $self->Bio::Graphics::Glyph::wiggle_data::_draw_coverage($summary,\@vals,@_);
+ $self->Bio::Graphics::Glyph::wiggle_data::_draw_coverage($summary,\@vals,@_);
}
}
}
@@ -240,6 +242,8 @@ sub draw_peaks {
my $flip = $self->{flip};
+ $self->{peak_cache} = []; # remember coordinates of the peaks
+
foreach my $peak (@peaks) {
my $x1 = $left + ($peak->{start} - $f_start) * $x_scale;
my $x2 = $left + ($peak->{stop} - $f_start) * $x_scale;
@@ -275,10 +279,12 @@ sub draw_peaks {
($x1,$x2) = ($x2,$x1);
}
- $self->filled_box($gd,int($x1+0.5),int($y1+0.5),int($x2+0.5),int($y2+0.5),$bgcolor,$bgcolor,0.5) if abs($y2-$y1) > 0;
+ my @rect = (int($x1+0.5),int($y1+0.5),int($x2+0.5),int($y2+0.5));
+ $self->filled_box($gd,@rect,$bgcolor,$bgcolor,0.5) if abs($y2-$y1) > 0;
$gd->setThickness($lw);
$gd->line(int($x1+0.5),int($y1+0.5),int($x2+0.5),int($y1+0.5),$color);
$gd->setThickness(1);
+ push @{$self->{peak_cache}},[$peak,@rect];
}
}
}
@@ -401,46 +407,16 @@ sub level { -1 }
# Need to override this so we have a nice image map for overlayed peaks
sub boxes {
- my $self = shift;
-
- return if $self->glyph_subtype eq 'density'; # No boxes for density plot
- my($left,$top,$parent) = @_;
-
- my $feature = $self->feature;
- my @result;
- my ($handle) = eval{$feature->get_tag_values('peak_type')};
-
- if (!$handle) {
- return wantarray ? () : \();
- }
-
- $parent ||=$self;
- $top += 0; $left += 0;
-
- if ($handle) {
- my @peaks = $self->peaks;
- $self->add_feature(@peaks);
-
- my $x_scale = $self->scale;
- my $panel_start = $self->panel->start;
- my $f_start = $feature->start > $panel_start
- ? $feature->start
- : $panel_start;
- for my $part ($self->parts) {
- my $x1 = $f_start < $part->{start} ? int(($part->{start} - $f_start) * $x_scale) : 1;
- my $x2 = int(($part->{stop} - $f_start) * $x_scale);
- my $y1 = 0;
- my $y2 = $part->height + $self->pad_top;
- $x2++ if $x1==$x2;
- next if $x1 <= 0;
- push @result,[$part->feature,
- $left + $x1,$top+$self->top+$self->pad_top+$y1,
- $left + $x2,$top+$self->top+$self->pad_top+$y2,
- $parent];
- }
- }
+ my $self = shift;
+ my($left,$top,$parent) = @_;
- return wantarray ? @result : \@result;
+ return if $self->glyph_subtype eq 'density'; # No boxes for density plot
+ my @boxes = $self->SUPER::boxes(@_);
+
+ if (my $rects = $self->{peak_cache}) {
+ push @boxes,[@$_,$parent] foreach @$rects;
+ }
+ return wantarray ? @boxes : \@boxes;
}
View
18 lib/Bio/Graphics/Glyph/wiggle_data.pm
@@ -3,7 +3,6 @@ package Bio::Graphics::Glyph::wiggle_data;
use strict;
use base qw(Bio::Graphics::Glyph::minmax);
use File::Spec;
-
sub minmax {
my $self = shift;
my $parts = shift;
@@ -157,7 +156,6 @@ sub datatype {
my $self = shift;
my $feature = $self->feature;
my ($tag,$value);
-
for my $t ('wigfile','wigdata','densefile','coverage') {
if (my ($v) = eval{$feature->get_tag_values($t)}) {
$value = $v;
@@ -178,7 +176,6 @@ sub get_parts {
my $feature = $self->feature;
my ($start,$end) = $self->effective_bounds($feature);
my ($datatype,$data) = $self->datatype;
-
return $self->subsample($data,$start,$end) if $datatype eq 'wigdata';
return $self->create_parts_from_wigfile($data,$start,$end) if $datatype eq 'wigfile';
return $self->create_parts_for_dense_feature($data,$start,$end) if $datatype eq 'densefile';
@@ -251,6 +248,9 @@ sub create_parts_from_summary {
sub create_parts_from_wigfile {
my $self = shift;
my ($path,$start,$end) = @_;
+ if (ref $path && $path->isa('Bio::Graphics::Wiggle')) {
+ return $self->create_parts_for_dense_feature($path,$start,$end);
+ }
$path = $self->rel2abs($path);
if ($path =~ /\.wi\w{1,3}$/) {
eval "require Bio::Graphics::Wiggle" unless Bio::Graphics::Wiggle->can('new');
@@ -325,7 +325,7 @@ sub draw_wigfile {
warn $@;
return $self->SUPER::draw(@_);
}
- $self->_draw_wigfile($feature,$wig,@_);
+ $self->_draw_wigfile($feature,$wigfile,@_);
}
sub draw_wigdata {
@@ -424,7 +424,15 @@ sub _draw_coverage {
sub _draw_wigfile {
my $self = shift;
my $feature = shift;
- my $wig = shift;
+ my $wigfile = shift;
+
+ $self->feature->remove_tag('wigfile') if $self->feature->has_tag('wigfile');
+ $self->feature->add_tag_value('wigfile',$wigfile);
+
+ eval "require Bio::Graphics::Wiggle" unless Bio::Graphics::Wiggle->can('new');
+ my $wig = ref $wigfile && $wigfile->isa('Bio::Graphics::Wiggle')
+ ? $wigfile
+ : eval { Bio::Graphics::Wiggle->new($wigfile) };
$wig->smoothing($self->get_smoothing);
$wig->window($self->smooth_window);
View
248 lib/Bio/Graphics/Glyph/wiggle_density.pm
@@ -1,5 +1,4 @@
package Bio::Graphics::Glyph::wiggle_density;
-
use strict;
use base qw(Bio::Graphics::Glyph::wiggle_data
Bio::Graphics::Glyph::box
@@ -85,6 +84,231 @@ sub draw {
}
}
+sub draw_coverage {
+ my $self = shift;
+ my $feature = shift;
+ my $array = shift;
+
+
+ if (! $array || ref($array) ne 'ARRAY'){
+ unshift(@_,$array);
+ my @arr = (eval{$feature->get_tag_values('coverage')});
+ $array = $arr[0];
+ } else {
+ $array = [split ',',$array] unless ref $array;
+ }
+ return unless @$array;
+
+ my ($gd,$left,$top) = @_;
+
+ my ($start,$end) = $self->effective_bounds($feature);
+ my $length = $end - $start + 1;
+ my $bases_per_bin = ($end-$start)/@$array;
+ my @parts;
+ my $samples = $length < $self->panel->width ? $length
+ : $self->panel->width;
+ my $samples_per_base = $samples/$length;
+
+ for (my $i=0;$i<$samples;$i++) {
+ my $offset = $i/$samples_per_base;
+ my $v = $array->[$offset/$bases_per_bin];
+ push @parts,$v;
+ }
+
+ my ($x1,$y1,$x2,$y2) = $self->bounds($left,$top);
+ $self->draw_segment($gd,
+ $start,$end,
+ \@parts,
+ $start,$end,
+ 1,1,
+ $x1,$y1,$x2,$y2);
+}
+
+sub draw_segment {
+ my $self = shift;
+ my ($gd,
+ $start,$end,
+ $seg_data,
+ $seg_start,$seg_end,
+ $step,$span,
+ $x1,$y1,$x2,$y2) = @_;
+
+ # clip, because wig files do no clipping
+ $seg_start = $start if $seg_start < $start;
+ $seg_end = $end if $seg_end > $end;
+
+ # figure out where we're going to start
+ my $scale = $self->scale; # pixels per base pair
+ my $pixels_per_span = $scale * $span + 1;
+ my $pixels_per_step = 1;
+ my $length = $end-$start+1;
+
+ # if the feature starts before the data starts, then we need to draw
+ # a line indicating missing data (this only happens if something went
+ # wrong upstream)
+ if ($seg_start > $start) {
+ my $terminus = $self->map_pt($seg_start);
+ $start = $seg_start;
+ $x1 = $terminus;
+ }
+ # if the data ends before the feature ends, then we need to draw
+ # a line indicating missing data (this only happens if something went
+ # wrong upstream)
+ if ($seg_end < $end) {
+ my $terminus = $self->map_pt($seg_end);
+ $end = $seg_end;
+ $x2 = $terminus;
+ }
+
+ return unless $start < $end;
+
+ # get data values across the area
+ my $samples = $length < $self->panel->width ? $length
+ : $self->panel->width;
+ my $data = ref $seg_data eq 'ARRAY' ? $seg_data
+ : $seg_data->values($start,$end,$samples);
+
+ # scale the glyph if the data end before the panel does
+ my $data_width = $end - $start;
+ my $data_width_ratio;
+ if ($data_width < $self->panel->length) {
+ $data_width_ratio = $data_width/$self->panel->length;
+ }
+ else {
+ $data_width_ratio = 1;
+ }
+ return unless $data && ref $data && @$data > 0;
+
+ my $min_value = $self->min_score;
+ my $max_value = $self->max_score;
+
+ my ($min,$max,$mean,$stdev) = $self->minmax($data);
+ unless (defined $min_value && defined $max_value) {
+ $min_value ||= $min;
+ $max_value ||= $max;
+ }
+
+ my $rescale = $self->option('autoscale') eq 'z_score';
+ my ($scaled_min,$scaled_max);
+ if ($rescale) {
+ my $bound = $self->z_score_bound;
+ $scaled_min = -$bound;
+ $scaled_max = +$bound;
+ } else {
+ ($scaled_min,$scaled_max) = ($min_value,$max_value);
+ }
+
+ my $t = 0; for (@$data) {$t+=$_}
+
+ # allocate colors
+ # There are two ways to do this. One is a scale from min to max. The other is a
+ # bipartite scale using one color range from zero to min, and another color range
+ # from 0 to max. The latter behavior is triggered when the config file contains
+ # entries for "pos_color" and "neg_color" and the data ranges from < 0 to > 0.
+
+ my $poscolor = $self->pos_color || $self->fgcolor;
+ my $negcolor = $self->neg_color || $self->bgcolor;
+
+ my $data_midpoint = $self->midpoint;
+ $data_midpoint = 0 if $rescale;
+ my $bicolor = $poscolor != $negcolor
+ && $scaled_min < $data_midpoint
+ && $scaled_max > $data_midpoint;
+
+ my ($rgb_pos,$rgb_neg,$rgb);
+ if ($bicolor) {
+ $rgb_pos = [$self->panel->rgb($poscolor)];
+ $rgb_neg = [$self->panel->rgb($negcolor)];
+ } else {
+ $rgb = $scaled_max > $scaled_min ? ([$self->panel->rgb($poscolor)] || [$self->panel->rgb($self->bgcolor)])
+ : ([$self->panel->rgb($negcolor)] || [$self->panel->rgb($self->bgcolor)]);
+ }
+
+
+ my %color_cache;
+
+ @$data = reverse @$data if $self->flip;
+
+ if (@$data <= $self->panel->width) { # data fits in width, so just draw it
+
+ $pixels_per_step = $scale * $step;
+ $pixels_per_step = 1 if $pixels_per_step < 1;
+ my $datapoints_per_base = @$data/$length;
+ my $pixels_per_datapoint = $self->panel->width/@$data * $data_width_ratio;
+
+ my %temps;
+ map{$temps{$_}++} (@$data);
+ my %colorss = ();
+ for (my $i = 0; $i <= @$data ; $i++) {
+ my $x = $x1 + $pixels_per_datapoint * $i;
+ my $data_point = $data->[$i];
+ defined $data_point || next;
+ $data_point = ($data_point-$mean)/$stdev if $rescale;
+ $data_point = $scaled_min if $scaled_min > $data_point;
+ $data_point = $scaled_max if $scaled_max < $data_point;
+ my ($r,$g,$b) = $bicolor
+ ? $data_point > $data_midpoint ? $self->calculate_color($data_point,$rgb_pos,
+ $data_midpoint,$scaled_max)
+ : $self->calculate_color($data_point,$rgb_neg,
+ $data_midpoint,$scaled_min)
+ : $self->calculate_color($data_point,$rgb,
+ $scaled_min,$scaled_max);
+
+ my $idx = $color_cache{$r,$g,$b} ||= $self->panel->translate_color($r,$g,$b);
+ $colorss{$idx} = $data_point;
+ $self->filled_box($gd,$x,$y1,$x+$pixels_per_datapoint,$y2,$idx,$idx);
+ }
+ (keys %colorss); # Alleviate a silent crash somewhere in GD that causes density graph get drawn as a solid-colored box
+ } else { # use Sheldon's code to subsample data
+ $pixels_per_step = $scale * $step;
+ my $pixels = 0;
+
+ # only draw boxes 2 pixels wide, so take the mean value
+ # for n data points that span a 2 pixel interval
+ my $binsize = 2/$pixels_per_step;
+ my $pixelstep = $pixels_per_step;
+ $pixels_per_step *= $binsize;
+ $pixels_per_step *= $data_width_ratio;
+ $pixels_per_span = 2;
+
+ my $scores = 0;
+ my $defined;
+
+ for (my $i = $start; $i < $end ; $i += $step) {
+ # draw the box if we have accumulated >= 2 pixel's worth of data.
+ if ($pixels >= 2) {
+ my $data_point = $defined ? $scores/$defined : 0;
+ $scores = 0;
+ $defined = 0;
+
+ $data_point = $scaled_min if $scaled_min > $data_point;
+ $data_point = $scaled_max if $scaled_max < $data_point;
+ my ($r,$g,$b) = $bicolor
+ ? $data_point > $data_midpoint ? $self->calculate_color($data_point,$rgb_pos,
+ $data_midpoint,$scaled_max)
+ : $self->calculate_color($data_point,$rgb_neg,
+ $data_midpoint,$scaled_min)
+ : $self->calculate_color($data_point,$rgb,
+ $scaled_min,$max_value);
+ my $idx = $color_cache{$r,$g,$b} ||= $self->panel->translate_color($r,$g,$b);
+ $self->filled_box($gd,$x1,$y1,$x1+$pixels_per_span,$y2,$idx,$idx);
+ $x1 += $pixels;
+ $pixels = 0;
+ }
+
+ my $val = shift @$data;
+ # don't include undef scores in the mean calculation
+ # $scores is the numerator; $defined is the denominator
+ $scores += $val if defined $val;
+ $defined++ if defined $val;
+
+ # keep incrementing until we exceed 2 pixels
+ # the step is a fraction of a pixel, not an integer
+ $pixels += $pixelstep;
+ }
+ }
+}
+
sub draw_plot {
my $self = shift;
my $parts = shift;
@@ -160,6 +384,28 @@ sub draw_plot {
return 1;
}
+sub _draw_coverage {
+ my $self = shift;
+ my $feature = shift;
+ my $array = shift;
+
+ $array = [split ',',$array] unless ref $array;
+ return unless @$array;
+
+ my ($start,$end) = $self->effective_bounds($feature);
+ my $bases_per_bin = ($end-$start)/@$array;
+ my $pixels_per_base = $self->scale;
+ my @parts;
+ for (my $pixel=0;$pixel<$self->width;$pixel++) {
+ my $offset = $pixel/$pixels_per_base;
+ my $s = $start + $offset;
+ my $e = $s+1; # fill in gaps
+ my $v = $array->[$offset/$bases_per_bin];
+ push @parts,[$s,$s,$v];
+ }
+ $self->Bio::Graphics::Glyph::wiggle_density::draw_plot(\@parts,@_);
+}
+
sub calculate_color {
my $self = shift;
my ($s,$rgb,$min_score,$max_score) = @_;
View
2 lib/Bio/Graphics/Glyph/wiggle_whiskers.pm
@@ -140,7 +140,7 @@ sub draw {
my $feature = $self->feature;
my $stats = eval {$feature->statistical_summary($self->width)};
- if ($@ =~ /can't locate object method/i) {
+ if ($@ =~ /can\'t locate object method/i) {
warn "This glyph only works properly with features that have a statistical_summary() method, but you passed a ",ref($feature)," object";
return;
}
View
84 lib/Bio/Graphics/Glyph/wiggle_xyplot.pm
@@ -98,6 +98,34 @@ sub draw {
return $result;
}
+sub draw_coverage {
+ my $self = shift;
+ my $feature = shift;
+ my $array = shift;
+ if (! $array || ref($array) ne 'ARRAY'){
+ unshift(@_,$array);
+ my @arr = (eval{$feature->get_tag_values('coverage')});
+ $array = $arr[0];
+ } else {
+ $array = [split ',',$array] unless ref $array;
+ }
+ return unless @$array;
+
+ my ($start,$end) = $self->effective_bounds($feature);
+ my $bases_per_bin = ($end-$start)/@$array;
+ my $pixels_per_base = $self->scale;
+ my @parts;
+ for (my $pixel=0;$pixel<$self->width;$pixel++) {
+ my $offset = $pixel/$pixels_per_base;
+ my $s = $start + $offset;
+ my $e = $s+1; # fill in gaps
+ my $v = $array->[$offset/$bases_per_bin];
+ push @parts,[$s,$s,$v];
+ }
+ $self->draw_plot(\@parts,@_);
+}
+
+
sub draw_plot {
my $self = shift;
my $parts = shift;
@@ -110,7 +138,7 @@ sub draw_plot {
# There is a minmax inherited from xyplot as well as wiggle_data, and I don't want to
# rely on Perl's multiple inheritance DFS to find the right one.
- my ($min_score,$max_score,$mean,$stdev) = $self->minmax($parts);
+ my ($min_score,$max_score,$mean,$stdev) = $self->minmax($parts);
my $rescale = $self->option('autoscale') eq 'z_score';
my $side = $self->_determine_side();
@@ -119,8 +147,8 @@ sub draw_plot {
$scaled_min = int(($min_score-$mean)/$stdev + 0.5);
$scaled_max = int(($max_score-$mean)/$stdev + 0.5);
my $bound = $self->z_score_bound;
- $scaled_max = $bound if $scaled_max > $bound;
- $scaled_min = -$bound if $scaled_min < -$bound;
+ $scaled_max = $bound if $scaled_max >= 0;
+ $scaled_min = -$bound if $scaled_min <= 0;
}
elsif ($side) {
$scaled_min = int($min_score - 0.5);
@@ -148,10 +176,12 @@ sub draw_plot {
my $y_origin = $scaled_min <= 0 && $pivot ne 'min' ? $bottom - (0 - $scaled_min) * $y_scale : $bottom;
$y_origin = int($y_origin+0.5);
+
$self->panel->startGroup($gd);
- $self->_draw_grid($gd,$x_scale,$scaled_min,$scaled_max,$dx,$dy,$y_origin) unless ($self->option('no_grid') == 1);
+ $self->_draw_grid($gd,$x_scale,$scaled_min,$scaled_max,$dx,$dy,$y_origin) unless $self->option('no_grid');
$self->panel->endGroup($gd);
+
return unless $scaled_max > $scaled_min;
my $lw = $self->linewidth;
@@ -286,6 +316,52 @@ sub draw_plot {
$self->panel->endGroup($gd);
}
+sub make_key_feature {
+ my $self = shift;
+ my $scale = 1/$self->scale; # base pairs/pixel
+
+ # one segments, at pixels 0->80
+ my $offset = $self->panel->offset;
+
+ my $start = 0 * $scale + $offset;
+ my $end = 80*$scale+$offset;
+
+ my $span = int(0.5+($end-$start)/50);
+
+ eval "require Bio::Graphics::Wiggle; 1";
+ my $wig = Bio::Graphics::Wiggle->new(undef,
+ 1,
+ {seqid=>'chr1',
+ start => int($start),
+ span => $span});
+
+ my @values = map {
+ (sin($_/60)+sin($_/12))*100+rand(100)
+ } 0..50;
+ my $min = $values[0];
+ my $max = $values[0];
+ my $tot = 0;
+ for (@values) {$min = $_ if $min > $_;
+ $max = $_ if $max < $_;
+ $tot += $_;
+ }
+ $wig->min($min);
+ $wig->max($max);
+ $wig->stdev(120); # just make it up
+ $wig->mean($tot/@values);
+
+ $wig->set_value(($_*$span)+$start=>$values[$_-1]) for (0..50);
+
+ my $feature =
+ Bio::Graphics::Feature->new(-start => $start,
+ -end => $end,
+ -name => $self->make_key_name(),
+ -strand => '+1',
+ -attributes => {wigfile=>$wig},
+ );
+ return $feature;
+}
+
1;
__END__
View
18 lib/Bio/Graphics/Panel.pm
@@ -1014,19 +1014,21 @@ sub _translate_color {
$colors[0] =~ /^rgb\((\d+),(\d+),(\d+)\)$/i
) {
my (@rgb) = map {/(\d+)%/ ? int(255 * $1 / 100) : $_} ($1,$2,$3);
- $index = $self->colorAllocateAlpha($gd,@rgb,$default_alpha);
+ $index = $gd->colorAllocateAlpha(@rgb,$default_alpha);
}
elsif ($colors[0] eq 'transparent') {
$index = $gd->colorAllocateAlpha(255,255,255,127);
}
elsif ($colors[0] =~ /^(\w+):([\d.]+)/) { # color:alpha
my @rgb = $self->color_name_to_rgb($1);
+ @rgb = (0,0,0) unless @rgb;
my $alpha = $self->adjust_alpha($2);
$index = $gd->colorAllocateAlpha(@rgb,$alpha);
}
elsif ($default_alpha < 127) {
my @rgb = $self->color_name_to_rgb($colors[0]);
- $index = $gd->colorAllocateAlpha(@rgb,$default_alpha);
+ @rgb = (0,0,0) unless @rgb;
+ $index = $gd->colorAllocateAlpha(@rgb,$default_alpha);
}
else {
$index = defined $table->{$colors[0]} ? $table->{$colors[0]} : 1;
@@ -1959,10 +1961,14 @@ object to the panel. This stylesheet will be called to determine both
the glyph and the glyph options. If both a stylesheet and direct
options are provided, the latter take precedence.
-The B<-color_series> argument, if true, causes the track to ignore
-the -bgcolor setting and instead to assign glyphs a series of
-contrasting colors. This is usually used in combination with
--bump=>'overlap' in order to create overlapping features.
+The B<-color_series> argument causes the track to ignore the -bgcolor
+setting and instead to assign glyphs a series of contrasting
+colors. This is usually used in combination with -bump=>'overlap' in
+order to create overlapping features. A true value activates the color
+series. You may adjust the default color series using the
+B<-color_cycle> option, which is either a reference to an array of
+Bio::Graphics color values, or a space-delimited string of color
+names/value.
If successful, add_track() returns an Bio::Graphics::Glyph object.
You can use this object to add additional features or to control the
View
16 lib/Bio/Graphics/Wiggle.pm
@@ -288,12 +288,7 @@ sub end {
return $self->{end};
}
-sub DESTROY {
- my $self = shift;
- if ($self->{dirty} && $self->{write}) {
- $self->_writeoptions($self->{options});
- }
-}
+sub DESTROY { shift->write }
sub erase {
my $self = shift;
@@ -344,6 +339,15 @@ sub smoothing {
$d;
}
+sub write {
+ my $self = shift;
+ if ($self->{dirty} && $self->{write}) {
+ $self->_writeoptions($self->{options});
+ undef $self->{dirty};
+ $self->fh->flush;
+ }
+}
+
sub _readoptions {
my $self = shift;
my $fh = $self->fh;
View
15 lib/Bio/Graphics/Wiggle/Loader.pm
@@ -8,6 +8,7 @@ package Bio::Graphics::Wiggle::Loader;
my $gff3_file = $loader->featurefile('gff3',$method,$source);
my $featurefile = $loader->featurefile('featurefile');
+ my @features = $loader->features();
=head1 USAGE
@@ -52,6 +53,11 @@ feature. In the case of "featurefile", the returned file will contain
GBrowse stanzas that describe a reasonable starting format to display
the data.
+=item @features = $loader->features
+
+Returns one or more Bio::Graphics::Features objects, which can be used to
+create Bio::Graphics tracks with the wiggle_xyplot (and related) glyphs.
+
=item $loader->allow_sampling(1)
If allow_sampling() is passed a true value, then very large files
@@ -109,6 +115,7 @@ use Statistics::Descriptive;
use IO::Seekable;
use File::Spec;
use Bio::Graphics::Wiggle;
+use Bio::Graphics::FeatureFile;
use Text::ParseWords();
use File::stat;
use CGI 'escape';
@@ -250,6 +257,7 @@ sub featurefile {
for my $seqid (@seqid) {
$seqid or next;
+ $tracks->{$track}{seqids}{$seqid}{wig}->write();
my $attributes = join ';',(@attributes,"wigfile=$seqids->{$seqid}{wigpath}");
if ($type eq 'gff3') {
push @lines,join "\t",($seqid,$source,$method,
@@ -271,6 +279,13 @@ sub featurefile {
return join("\n",@lines)."\n";
}
+sub features {
+ my $self = shift;
+ my $text = $self->featurefile('featurefile');
+ my $file = Bio::Graphics::FeatureFile->new(-text=>$text);
+ return $file->features;
+}
+
sub load {
my $self = shift;
View
14 scripts/glyph_help.pl
@@ -10,6 +10,7 @@
use Bio::Graphics::Wiggle;
my $MANUAL = 0;
+my $POD = 0;
my $LIST = 0;
my $PICT = 0;
my $VIEW = 0;
@@ -24,6 +25,8 @@
Options:
-m --manual Print the full manual page for the glyph, followed
by a summary of its options.
+ -r --raw Print the quick summary of the glyph\'s options in raw POD
+ format.
-l --list List all glyphs that are available for use.
-p --picture Create a PNG picture of what the indicated glyph looks like.
The PNG will be written to stdout
@@ -46,6 +49,7 @@
USAGE
GetOptions ('manual' => \$MANUAL,
+ 'raw' => \$POD,
'list' => \$LIST,
'picture' => \$PICT,
'view' => \$VIEW,
@@ -76,9 +80,13 @@
exit 0;
}
-system "perldoc",$class if $MANUAL;
-$class->options_man();
-
+if ($MANUAL) {
+ system "perldoc",$class;
+} elsif ($POD) {
+ print $class->options_pod();
+} else {
+ print $class->options_man();
+}
exit 0;
sub print_list {
View
384 scripts/index_cov_files.pl
@@ -0,0 +1,384 @@
+#!/usr/bin/env perl
+#
+# index_cov_files.pl
+#
+# [ to see proper formatting set tab==2 ]
+#
+# 2009-2010 Victor Strelets, FlyBase.org
+
+##$testing= 300000;
+
+#$debug= 1; $do_only_chr= '2L';
+#$debug= 1; $do_only_chr= '2L'; $do_only_subset= 'BS40_all_unique';
+#$debug= 1; $do_only_subset= 'Female_Heads'; $do_only_chr= '2';
+
+$apply_log= 0;
+$log2= log(2);
+%LogTabs= ( 0, 0, 1, 0 ); # for our purposes, these start values are right
+$log_magnifier= 1.0;
+
+ if( @ARGV && $ARGV[0]=~/^\-?log/i ) { $apply_log= 1; print "applying log\n"; }
+
+ print "\nIndexing COV files for use with the fb_shmiggle glyph\n";
+
+ eval {use lib("/LS/common/system-local/perl/lib"); };
+ use BerkeleyDB;
+
+ my $filesmask= '*.cov*';
+ $filemask= shift(@ARGV) if @ARGV;
+ $filemask=~s/\./\\\./g;
+ $filemask=~s/[\*]/\\*/g;
+ indexfeatdir('./',$filesmask) ;
+
+ exit();
+
+#*************************************************************************
+#
+#*************************************************************************
+
+sub indexfeatdir
+{
+ my($dir,$mask)= @_;
+
+ local(*D); opendir(D, $dir) || warn "can't open $dir";
+ my @files= grep( /\.(cov|wig)/i, readdir(D));
+ closedir(D);
+
+ my $datfilename= 'data.cat';
+ system("rm $datfilename") if -e $datfilename;
+ open(OUTDATF,'>'.$datfilename) || die "Cannot open $datfilename!";
+ my $bdbfilename= 'index.bdbhash';
+ unlink($bdbfilename) if( -e $bdbfilename );
+ %ResIndexHash= (); # !!GLOBAL
+ tie %ResIndexHash, "BerkeleyDB::Hash", -Filename => $bdbfilename, -Flags => DB_CREATE;
+
+ $max_signal = 0;
+ @SubsetNames= ();
+ foreach my $file (sort @files) {
+ indexCoverageFile($file); # coverage files are in fact wiggle files..
+ }
+ $ResIndexHash{'subsets'}= join("\t",@SubsetNames); # record subsets, just in case..
+ $ResIndexHash{'max_signal'}= $max_signal;
+ my @all_keys= keys %ResIndexHash;
+ foreach my $kkey ( sort @all_keys ) { print "\t$kkey => ".$ResIndexHash{$kkey}."\n"; }
+ if( $max_signal>10000 ) { print "WARNING: max_signal=$max_signal - TOO HIGH!! Re-run with '-log' option\n"; }
+ untie %ResIndexHash;
+ chmod(0666,$bdbfilename); # ! sometimes very important
+ close(OUTDATF);
+
+ return;
+}
+
+#*************************************************************************
+#
+#*************************************************************************
+
+sub indexCoverageFile
+{
+ my($file)= @_;
+ if( $debug && defined $do_only_subset ) {
+ return unless $file=~/^${do_only_subset}\./; }
+ my $zcat= get_zcat($file);
+ local(*INF);
+ open(INF,"$zcat $file |") || die "Can't open $file";
+ print "\t$file\n";
+ my $SubsetName= ($file=~/^([^\.]+)\./) ? $1 : $file;
+ push(@SubsetNames,$SubsetName);
+ my $chromosome= "";
+ my @offsets= ();
+ # following setting is very important for performance (in some cases)
+ # value 1000 (otherwise good) on K.White dataset was causing start of reading 100K before the actually required point..
+ my $step= 1000; # step in coverage file lines ()signal reads to save start-offset
+ my $coordstep= 20000; # step in coords to save start-offset
+ my $counter= 0;
+ my $offset= tell(OUTDATF);
+ $ResIndexHash{$SubsetName}= $offset; # record offset where new subset data starts
+ my $old_signal= 0;
+ my $oldcoord= -200000;
+ my $FileFormat= 1;
+ my $StartCoord= 0;
+ my $lastRecordedCoord= -200000;
+ while( (my $str= <INF>) ) {
+ $offset= tell(OUTDATF);
+
+ # correct variant of GEO preferred subset spec
+ if( $str=~m/^(track[ \t]+type=wiggle_0)\s*\n$/i ) { # new subset starting
+ $str= $1 . ' name="' . $SubsetName ."\"\n"; }
+
+ # following is a GEO preferred subset spec
+ if( $str=~m/^track[ \t]+type=wiggle_0[ \t]+name="([^"]+)"/i ) { # new subset starting
+ $FileFormat= 4;
+ $SubsetName= $1;
+ #$chromosome= "";
+ next; # because it is not a signal, should not be printed in this data loop
+ }
+
+ # fix for K.White files
+ elsif( $str=~m/^track[ \t]+type=bedGraph[ \t]+name="([^"]+)"/i ) { # new subset starting
+ $FileFormat= 4;
+ $SubsetName=~s/_(combined|coverage)$//i; $SubsetName=~s/_(combined|coverage)$//i; $SubsetName=~s/^G[A-Z]{2}\d+[_\-]//;
+ #$chromosome= "";
+ next; # because it is not a signal, should not be printed in this data loop
+ }
+
+ elsif( $str=~m/^variableStep[ \t]+(chr(om(osome)?)?|arm)=(\w+)/i ) { # potentially new arm starting
+ $FileFormat= 4;
+ my $new_chromosome= $4;
+ $new_chromosome=~s/^chr(omosome)?//i;
+ if( $new_chromosome ne $chromosome ) {
+ dumpOffsets($SubsetName.':'.$chromosome,@offsets) unless $chromosome eq ""; # previous subset:arm
+ $chromosome= $new_chromosome;
+ print OUTDATF "# subset=$SubsetName chromosome=$chromosome\n";
+ $offset= tell(OUTDATF);
+ $ResIndexHash{$SubsetName.':'.$chromosome}= $offset; # record offset where new subset:arm data starts
+ @offsets= ("-200000\t$offset");
+ print OUTDATF "-200000\t0\n"; # insert one fictive zero read
+ $offset= tell(OUTDATF);
+ print OUTDATF "0\t0\n"; # insert one more fictive zero read
+ push(@offsets,"0\t$offset");
+ print "\t\t$SubsetName:$chromosome\n";
+ $counter= 0; $old_signal= 0; $oldcoord= 0; $lastRecordedCoord= 0;
+ }
+ next; # because it is not a signal, should not be printed in this data loop
+ }
+
+ elsif( $str=~m/^FixedStep[ \t]+(chr(om(osome)?)?|arm)=(\w+)[ \t]+Start=(\d+)/i ) { # potentially new arm starting
+ $FileFormat= 3;
+ my $new_chromosome= $4;
+ $StartCoord= $5;
+ $new_chromosome=~s/^chr(omosome)?//i;
+ if( $new_chromosome ne $chromosome ) {
+ dumpOffsets($SubsetName.':'.$chromosome,@offsets) unless $chromosome eq ""; # previous subset:arm