Skip to content

HTTPS clone URL

Subversion checkout URL

You can clone with
or
.
Download ZIP
branch: master
Fetching contributors…

Cannot retrieve contributors at this time

executable file 256 lines (178 sloc) 5.764 kB
#!/usr/local/bin/perl
=head1 NAME
gff2ps - you will want to change this script
=head2 SYNOPSIS
perl gff2ps < file.gff > file.ps
=head2 DESCRIPTION
This script provides GFF to postscript handling. Due to the ... ummm
... potential for flexible reinterpretation that is GFF, this script
will almost certainly need modifying for anyone elses use (basically,
you need to know what you want to get out of the GFF file and how to
draw it). But it does include code to draw the most challenging thing
out there - genes - and should give you a good example of where to
start
=head2 AUTHOR
Ewan Birney
=cut
use Bio::Tools::GFF;
my $font = 8;
my $scale = 200;
my $rotate = 1;
my $feature_off = 0;
use Getopt::Long;
&GetOptions(
"scale=i" => \$scale,
"font=i" => \$font,
"rotate=i" => \$rotate,
"start=i" => \$feature_off
);
my $gffio = Bio::Tools::GFF->new(-fh => \*STDIN, -gff_version => 1);
my $feature;
use Data::Dumper;
my %set;
# loop over the input stream
while( my $f = $gffio->next_feature()) {
$f->start($f->start - $feature_off);
$f->end ($f->end - $feature_off);
if( $f->start < 0 ) {
next;
}
if( $f->start > $scale*1000 ) {
next;
}
if( $f->primary_tag =~ /coding_exon/ ) {
#print STDERR "Seen ",$f->start," ",$f->end,"\n";
($group) = $f->each_tag_value('group');
$group =~ s/\s+//g;
if( !defined $set{$group} ) {
$set{$group} = Bio::SeqFeature::Generic->new();
$set{$group}->seqname($f->seqname);
$set{$group}->primary_tag('transcript');
$set{$group}->source_tag($f->source_tag);
$set{$group}->add_tag_value('id',$group);
}
$set{$group}->add_sub_SeqFeature($f,'EXPAND');
$set{$group}->strand($f->strand);
}
}
$gffio->close();
#foreach my $set ( values %set ) {
# print $set->gff_string,"\n";
# foreach $sub ( $set->sub_SeqFeature ) {
# print $sub->gff_string,"\n";
# }
#}
# sort into forward and reverse strands
my @forward;
my @reverse;
$max = 0;
foreach my $set ( values %set ) {
if( $set->end > $max ) {
$max = $set->end;
}
if( $set->strand == -1 ) {
push(@reverse,$set);
} else {
push(@forward,$set);
}
}
@forward = sort { $a->start <=> $b->start } @forward;
@reverse = sort { $a->start <=> $b->start } @reverse;
&print_header(\*STDOUT);
if( $rotate ) {
print "0 700 translate\n";
print "-90 rotate\n";
}
print "0 200 moveto 900 200 lineto stroke\n";
my $bp_max = $scale*900;
for(my $bp = 0;$bp < $bp_max ;$bp = $bp + 5000) {
print STDOUT $bp/$scale," 200 moveto ",$bp/$scale," 197 lineto\n";
$text = int( $feature_off + ($bp/1000));
print STDOUT $bp/$scale," 195 moveto ($text) show\n";
}
&draw_gene(\@forward,1,$scale,220,\*STDOUT);
&draw_gene(\@reverse,-1,$scale,180,\*STDOUT);
print "showpage\n";
sub draw_gene {
my ($gene_array,$strand,$scale,$offset,$fh) = @_;
my @bump_array;
my $bump_row_max = 1;
my $bump_end = int $max/$scale;
$bump_array[0] = '0' x $bump_end;
foreach my $f ( @$gene_array ) {
#
# Bump it baby!
#
# We keep an array of strings for currently draw areas. Do this in pixel
# coordinates to save on space. If the region has all 0's then we know we
# can draw here. If so, we set it to 1's. If not, we go up a row and see if
# we can fit it in there. If we exhausted the rows we make a new row.
$bump_start = (int $f->start/$scale)-1;
$bump_len = int(($f->end - $f->start)/$scale) +1;
# text might be longer than gene. Mystic number 5 looks good for 8 point helvetica
# you will have to change this otherwise.
my ($gene_id) = $f->each_tag_value('id');
if( (length $gene_id)*5 > $bump_len ) {
$bump_len = (length $gene_id)*5;
}
# figure out the first place to fit in this gene;
for($i=0;$i<$bump_row_max;$i++) {
#print STDERR "Seeing $bump_start $bump_len $i ",substr($bump_array[$i],$bump_start,$bump_len),"\n";
if( substr($bump_array[$i],$bump_start,$bump_len) !~ /1/ ) {
#print STDERR "Going to break with $i\n";
last;
}
}
#print STDERR "i is $i\n";
# if $i == bump_row_max then we need a new bump row
if( $i == $bump_row_max ) {
$bump_array[$bump_row_max] = '0' x $bump_end;
$bump_row_max++;
}
# now blank out this bump row to 1's
substr($bump_array[$i],$bump_start,$bump_len) = '1' x $bump_len;
# now print it out ;)
#
# Need to be portable between strands. Gene hats go the
# other way up on reverse strand, but not the text.
#
if( $strand == 1 ) {
$text = $offset+($i*20)+1;
$bottom = $offset+($i*20)+10;
$top = $offset+($i*20)+20;
$mid = $offset+($i*20)+15;
} else {
$text = $offset-($i*20)-19;
$bottom = $offset-($i*20);
$top = $offset-($i*20)-10;
$mid = $offset-($i*20)-5;
}
print $fh $f->start/$scale," ",$text," moveto\n";
print $fh "($gene_id) show\n";
my $prev = undef;
foreach $exon ( $f->sub_SeqFeature ) {
print $fh $exon->start/$scale," ",$bottom," moveto\n";
print $fh $exon->end/$scale," ",$bottom," lineto\n";
print $fh $exon->end/$scale," ",$top, " lineto\n";
print $fh $exon->start/$scale," ",$top," lineto\n";
print $fh "closepath stroke\n";
# draw the intron hat
if( defined $prev ) {
print $prev->end/$scale," ",$mid," moveto\n";
my $intron_len = $exon->start - $prev->end;
print $fh ($prev->end+($intron_len/2))/$scale," ",$top," lineto\n";
print $fh $exon->start/$scale," ",$mid," lineto stroke\n";
}
$prev = $exon;
}
}
}
sub print_header {
my $fh = shift;
print $fh <<EOF;
%!PS-Adobe-2.0
% Created by Genome2ps. Ewan Birney <birney\@ebi.ac.uk>
0.5 setlinewidth
/Helvetica findfont $font scalefont setfont
EOF
}
Jump to Line
Something went wrong with that request. Please try again.