examples/tools/gff2ps.pl

   1 #!/usr/local/bin/perl
   2
   3
   4 =head1 NAME
   5
   6 gff2ps - you will want to change this script
   7
   8 =head2 SYNOPSIS
   9
  10    perl gff2ps < file.gff > file.ps
  11
  12 =head2 DESCRIPTION
  13
  14 This script provides GFF to postscript handling. Due to the ... ummm
  15 ... potential for flexible reinterpretation that is GFF, this script
  16 will almost certainly need modifying for anyone elses use (basically,
  17 you need to know what you want to get out of the GFF file and how to
  18 draw it). But it does include code to draw the most challenging thing
  19 out there - genes - and should give you a good example of where to
  20 start
  21
  22 =head2 AUTHOR
  23
  24 Ewan Birney
  25
  26 =cut
  27
  28
  29 use Bio::Tools::GFF;
  30
  31 my $font      = 8;
  32 my $scale     = 200;
  33 my $rotate    = 1;
  34 my $feature_off = 0;
  35
  36 use Getopt::Long;
  37
  38 &GetOptions(
  39             "scale=i"   => \$scale,
  40             "font=i"    => \$font,
  41             "rotate=i"  => \$rotate,
  42             "start=i"   => \$feature_off
  43             );
  44
  45
  46 my $gffio = Bio::Tools::GFF->new(-fh => \*STDIN, -gff_version => 1);
  47 my $feature;
  48
  49 use Data::Dumper;
  50
  51 my %set;
  52
  53 # loop over the input stream
  54 while( my $f = $gffio->next_feature()) {
  55     $f->start($f->start - $feature_off);
  56     $f->end  ($f->end   - $feature_off);
  57
  58     if( $f->start < 0 ) {
  59         next;
  60     }
  61
  62     if( $f->start > $scale*1000 ) {
  63       next;
  64     }
  65
  66
  67     if( $f->primary_tag =~ /coding_exon/ ) {
  68         #print STDERR "Seen ",$f->start," ",$f->end,"\n";
  69         ($group) = $f->each_tag_value('group');
  70         $group =~ s/\s+//g;
  71
  72
  73
  74         if( !defined $set{$group} ) {
  75             $set{$group} = Bio::SeqFeature::Generic->new();
  76             $set{$group}->seqname($f->seqname);
  77             $set{$group}->primary_tag('transcript');
  78             $set{$group}->source_tag($f->source_tag);
  79             $set{$group}->add_tag_value('id',$group);
  80         }
  81         $set{$group}->add_sub_SeqFeature($f,'EXPAND');
  82         $set{$group}->strand($f->strand);
  83     }
  84 }
  85 $gffio->close();
  86
  87
  88 #foreach my $set ( values %set ) {
  89 #    print $set->gff_string,"\n";
  90 #    foreach $sub ( $set->sub_SeqFeature ) {
  91 #       print $sub->gff_string,"\n";
  92 #    }
  93 #}
  94
  95
  96 # sort into forward and reverse strands
  97
  98 my @forward;
  99 my @reverse;
 100
 101 $max = 0;
 102
 103 foreach my $set ( values %set ) {
 104     if( $set->end > $max ) {
 105         $max = $set->end;
 106     }
 107
 108     if( $set->strand == -1 ) {
 109         push(@reverse,$set);
 110     } else {
 111         push(@forward,$set);
 112     }
 113 }
 114
 115 @forward = sort { $a->start <=> $b->start } @forward;
 116 @reverse = sort { $a->start <=> $b->start } @reverse;
 117
 118 &print_header(\*STDOUT);
 119
 120 if( $rotate ) {
 121    print "0 700 translate\n";
 122    print "-90 rotate\n";
 123 }
 124
 125 print "0 200 moveto 900 200 lineto stroke\n";
 126
 127 my $bp_max = $scale*900;
 128
 129 for(my $bp = 0;$bp < $bp_max ;$bp = $bp + 5000) {
 130     print STDOUT $bp/$scale," 200 moveto ",$bp/$scale," 197 lineto\n";
 131     $text = int( $feature_off + ($bp/1000));
 132     print STDOUT $bp/$scale," 195 moveto ($text) show\n";
 133 }
 134
 135
 136 &draw_gene(\@forward,1,$scale,220,\*STDOUT);
 137 &draw_gene(\@reverse,-1,$scale,180,\*STDOUT);
 138
 139 print "showpage\n";
 140
 141
 142
 143
 144 sub draw_gene {
 145     my ($gene_array,$strand,$scale,$offset,$fh) = @_;
 146
 147
 148     my @bump_array;
 149     my $bump_row_max = 1;
 150     my $bump_end = int $max/$scale;
 151
 152     $bump_array[0] = '0' x $bump_end;
 153
 154
 155     foreach my $f ( @$gene_array ) {
 156
 157         #
 158         # Bump it baby!
 159         #
 160
 161         # We keep an array of strings for currently draw areas. Do this in pixel
 162         # coordinates to save on space. If the region has all 0's then we know we
 163         # can draw here. If so, we set it to 1's. If not, we go up a row and see if
 164         # we can fit it in there. If we exhausted the rows we make a new row.
 165
 166         $bump_start = (int $f->start/$scale)-1;
 167         $bump_len   = int(($f->end - $f->start)/$scale) +1;
 168
 169         # text might be longer than gene. Mystic number 5 looks good for 8 point helvetica
 170         # you will have to change this otherwise.
 171
 172         my ($gene_id) = $f->each_tag_value('id');
 173         if( (length $gene_id)*5 > $bump_len ) {
 174             $bump_len = (length $gene_id)*5;
 175         }
 176
 177         # figure out the first place to fit in this gene;
 178         for($i=0;$i<$bump_row_max;$i++) {
 179             #print STDERR "Seeing $bump_start $bump_len $i ",substr($bump_array[$i],$bump_start,$bump_len),"\n";
 180
 181             if( substr($bump_array[$i],$bump_start,$bump_len) !~ /1/ ) {
 182                 #print STDERR "Going to break with $i\n";
 183                 last;
 184             }
 185         }
 186         #print STDERR "i is $i\n";
 187         # if $i == bump_row_max then we need a new bump row
 188         if( $i == $bump_row_max ) {
 189             $bump_array[$bump_row_max] = '0' x $bump_end;
 190             $bump_row_max++;
 191         }
 192
 193         # now blank out this bump row to 1's
 194
 195         substr($bump_array[$i],$bump_start,$bump_len) = '1' x $bump_len;
 196
 197         # now print it out ;)
 198
 199
 200         #
 201         # Need to be portable between strands. Gene hats go the
 202         # other way up on reverse strand, but not the text.
 203         #
 204
 205         if( $strand == 1 ) {
 206             $text   = $offset+($i*20)+1;
 207             $bottom = $offset+($i*20)+10;
 208             $top    = $offset+($i*20)+20;
 209             $mid    = $offset+($i*20)+15;
 210         } else {
 211             $text   = $offset-($i*20)-19;
 212             $bottom = $offset-($i*20);
 213             $top    = $offset-($i*20)-10;
 214             $mid    = $offset-($i*20)-5;
 215         }
 216
 217         print $fh $f->start/$scale," ",$text," moveto\n";
 218         print $fh "($gene_id) show\n";
 219
 220         my $prev = undef;
 221
 222         foreach $exon ( $f->sub_SeqFeature ) {
 223             print $fh $exon->start/$scale," ",$bottom," moveto\n";
 224             print $fh $exon->end/$scale," ",$bottom," lineto\n";
 225             print $fh $exon->end/$scale," ",$top, " lineto\n";
 226             print $fh $exon->start/$scale," ",$top," lineto\n";
 227             print $fh "closepath stroke\n";
 228
 229             # draw the intron hat
 230             if( defined $prev ) {
 231                 print $prev->end/$scale," ",$mid," moveto\n";
 232                 my $intron_len = $exon->start - $prev->end;
 233
 234                 print $fh ($prev->end+($intron_len/2))/$scale," ",$top," lineto\n";
 235                 print $fh $exon->start/$scale," ",$mid," lineto stroke\n";
 236             }
 237
 238             $prev = $exon;
 239         }
 240
 241     }
 242 }
 243
 244
 245 sub print_header {
 246     my $fh = shift;
 247
 248     print $fh <<EOF;
 249 %!PS-Adobe-2.0
 250 % Created by Genome2ps. Ewan Birney <birney\@ebi.ac.uk>
 251 0.5 setlinewidth
 252 /Helvetica findfont $font scalefont setfont
 253 EOF
 254
 255 }