#!/usr/local/bin/perl5.6.0 -w

# $Header: /nfs/disk222/yeastpub/Repository/zoo/general/embl_to_ps.pl,v 1.30 2000/10/05 13:14:31 kmr Exp $

# purpose: This script will generate an overview of embl entry in postscript.
# purpose: See the usage message for details on the arguments to the script. 

# bugs:    "Axiom:        Any non-trivial program contains some bugs
# bugs:     Corollary:    If a program has no bugs, it is useless"

use strict;

use AnT::Analyser;
use AnT::Feature;
use AnT::Range;
use Diagram::Data;
use Pathogens::Artemis;

# the number of bases on each line of the output
my $line_base_length = 20000;

# an page is 72*8.25 = 594 wide and 72*11 = 792 high

my $page_width = 72 * 8.25;
my $border_width = 72 * 0.5;
# this is the length of the line in points
my $page_height = 72*11;
my $page_height_less_border = $page_height - $border_width * 2;
my $line_spacing = $page_height_less_border / 15;
my $feature_height = $page_height_less_border / 200;
my $scale_match_distance = $feature_height;
my $scale_feature_distance = $feature_height;
my $label_feature_distance = $feature_height / 3;
my $font_size = 4.000;
my $tick_mark_spacing = 1000;  # bases
my $thick_mark_spacing = 10000;  # bases
my $label_angle = 45;
my $scale_number_distance = 3;
my $scale_font_size = 6.0;
my $shade_spacing = 2;
my $shade_right_keys = "";
my $shade_left_keys = "";
my $protein_features = "CDS,tRNA,rRNA,misc_RNA";
my $dna_line_features = "repeat_region,repeat_unit";
my $base_label_position = "both";
my $horizontal_labels = "false";
my $dna_feature_labels = "on_dna";
my $protein_feature_labels = "yes";
my $filter_in = "";
my $thick_tick_line_height = $feature_height;
my $tick_line_height = $feature_height * 0.75;

sub usage
{
  die <<EOF;
usage $0 -h (for help) 
   or $0 embl_file
   or $0 [-option=value ...] embl_file
options:
  (note that one unit of height or width is 1/72 of an inch)
  label_angle            - the angle of the feature label (0 is horizontal,
                           default is 45)
  feature_height         - the height of the each feature (default: 3.6)
  line_base_length       - the number of bases of each line (default: 20000)
  page_width             - the width of the page (default is 8.25 inches * 72 =
                           594)
  border_width           - the amount of space to leave at the edge of the page
                           (default is 36 = 0.5 inches)
  page_height            - the height of the page (default: 11 inches * 72 = 
                           792)
  line_spacing           - the distance between the dna lines (default: 48 = 
                           2/3 inches)
  scale_feature_distance - the distance between the scale (dna) line and the
                           feature (default: 3.6)
  font_size              - the amount by which the font should be scaled
                           (default: 4.0)
  scale_font_size        - the amount by which the font used to label the
                           dna scale line should be scaled (default: 8.0)
  scale_number_distance  - the distance between the dna scale line and the
                           labels at each end (default: 3)
  label_feature_distance - the distance between the feature and it's label
                           feature (default: 3.6)
  tick_mark_spacing      - the distance (in bases) between the small tick marks
                           on the scale line (default: 1000)
  thick_mark_spacing     - the distance (in bases) between the large tick marks
                           on the scale line (default: 10000)
  shade_spacing          - the distance between the diagonal lines used for
                           shading
  shade_left_keys        - the keys of the features that should be shaded with
                           left-slanting lines.  (comma separated list)
  shade_right_keys       - the keys of the features that should be shaded with
                           right-slanting lines.  Keys can be in this list and
                           the shade_left_keys list.  (comma separated list)
  protein_features       - the keys of the features that should be drawn beside
                           the dna line.  (comma separated list) (default:
                           CDS,tRNA,rRNA,misc_RNA)
  dna_line_features      - the keys of the features that should be drawn on the
                           dna line.  (comma separated list) (default:
                           repeat_region,repeat_unit)
  base_label_position    - this controls which end of the dna line should have
                           the base count label.  Possible values are none, 
                           left right and both.  (default: both) 
  horizontal_labels      - if true the labels on the protein features will be
                           horizontal and centred, otherwise the labels will
                           be drawn at the angle given by label_angle.
                           (default: false) 
  dna_feature_labels     - if this option is set to "on_dna", DNA line features
                           will have labels centred on the DNA line.  If this 
                           option is "no" no labels will be drawn.  Otherwise
                           the label will be drawn above the feature.
                           (default: on_dna)
  protein_feature_labels - if this option is set to "no" then no labels will be
                           drawn on protein features, otherwise the label
                           will be drawn.
                           (default: yes)

Note: Those features that have /note="signal" as a qualifier will have
a marker drawn near the start of the first exon.

EOF

}

if (@ARGV == 0 || $ARGV[0] eq "-h") {
  usage ();
}

while ($ARGV[0] =~ /^-(.*)=(.*)/) {
  shift;
  my $var_name = $1;
  my $var_value = quotemeta $2;
  
  # set $font_size if the option is -font_size=something

  eval <<EOF;
no strict;
if (defined \$$var_name) {
  \$$var_name = "$var_value"
} else {
  die "$0: unknown option: $var_name\n";
}
EOF

  if (defined $@ && $@) {
    print STDERR "$@\n";
    usage ();
  } else {
    next;
  }
}

if ($ARGV[0] =~ /^-(.*)/) {
  print STDERR "$0: unknown option: -$1\n";
  usage ();
}

if (@ARGV != 1) {
  usage ();
}

if ($horizontal_labels eq "false") {
  $horizontal_labels = 0;
}

my $line_length = $page_width - $border_width * 2;

sub get_features_in_range
{
  my $entry = shift;
  my $range = shift;

  my @return_vector = ();

  my @features = $entry->features ();

# loop over the features in the current entry
  for my $f (@features) {

    my @ranges = $f->ranges ();
    
    for (my $range_index = 0 ; $range_index < @ranges ; ++$range_index) {
      my $this_range = $ranges[$range_index];

      if ($this_range->overlaps ($range)) {
        push @return_vector, $f;
        last;
      }
    }
  }
  
  return @return_vector;
}


sub draw_exon
{
  my ($complement_flag, $dna_line_feature,
      $scale_line_y_pos, $line_start_base, $exon, $right_exon, $exon_index,
      $f) = @_;

  my $identifier = join ",", map {$_->start . ".." . $_->end} $f->ranges;
    
  my $x_start = 
    ($exon->start - $line_start_base) / $line_base_length * $line_length +
    $border_width;
  my $x_end = 
    ($exon->end - $line_start_base) / $line_base_length * $line_length +
    $border_width;
  my $y_top;
  my $y_bottom;

  my $match_y_top = 0;
  my $match_y_bottom = 0;

  if ($dna_line_feature) {
    $y_bottom = $scale_line_y_pos - $feature_height / 2;
    $y_top = $scale_line_y_pos + $feature_height / 2;
  } else {
    if ($complement_flag) {
      $y_bottom =
        $scale_line_y_pos - $feature_height - $scale_feature_distance;
      $y_top = $scale_line_y_pos - $scale_feature_distance;
      $match_y_bottom =
        $scale_line_y_pos - $feature_height / 2 - $scale_match_distance;
      $match_y_top = $scale_line_y_pos - $scale_match_distance;
    } else {
      $y_top = 
        $scale_line_y_pos + $feature_height + $scale_feature_distance;
      $y_bottom = $scale_line_y_pos + $scale_feature_distance;
      $match_y_top = 
        $scale_line_y_pos + $feature_height / 2 + $scale_match_distance;
      $match_y_bottom = $scale_line_y_pos + $scale_match_distance;
    }
  }



  print get_feature_colour ($f), " setrgbcolor\n";

  print <<EOF;

% feature $identifier
0.25 setlinewidth
newpath
$x_start $y_top M
$x_start $y_bottom L
$x_end $y_bottom L
$x_end $y_top L
closepath
fill
EOF
  
  # draw a black border
  print <<EOF;
0 0 0 setrgbcolor
newpath
$x_start $y_top M
$x_start $y_bottom L
$x_end $y_bottom L
$x_end $y_top L
closepath
stroke
EOF

  my $hp_match = $f->hp_match;

  if (defined $hp_match && $hp_match ne '"-"') {
    if ($hp_match =~ /".*\s+([\d\.]+)\%\s*(\(.*\))?\s*"/) {

      print 1, " ",(100-$1)/100.0, " ", (100-$1)/100.0, " setrgbcolor\n";

      print <<EOF;

% feature $identifier
0.25 setlinewidth
newpath
$x_start $match_y_top M
$x_start $match_y_bottom L
$x_end $match_y_bottom L
$x_end $match_y_top L
closepath
fill
EOF
      
    } else {
      die "can't understand this /hp_match qualifier: /hp_match=$hp_match\n";
    }
  }

  my $ml_unique = $f->has_q("ml_unique");

  print STDERR "$ml_unique\n";


  if ($ml_unique) {
    print <<EOF;
% feature $identifier

1 0 0 setrgbcolor
0.25 setlinewidth
newpath
$x_start $match_y_top M
$x_start $match_y_bottom L
$x_end $match_y_bottom L
$x_end $match_y_top L
closepath
fill
EOF
  }

  my $note = $f->note;

  if (defined $note && $note eq "\"signal\"" && $exon_index == 0) {
    # the height of each part of the marker
    my $line_height = $feature_height / 4;

    # draw a marker at the start of the exon
    
    print <<EOF;
0.4 setlinewidth
newpath
$x_start $y_top M
EOF

    if ($complement_flag) {
      print <<EOF;
-$line_height $line_height rmoveto
0 $line_height rmoveto
EOF
      for (my $i = 0 ; $i < 4 ; ++$i) {
        print <<EOF;
-$line_height -$line_height rlineto
$line_height -$line_height rlineto
EOF
      }
    } else {
      print <<EOF;
$line_height $line_height rmoveto
0 $line_height rmoveto
EOF
      for (my $i = 0 ; $i < 4 ; ++$i) {
        print <<EOF;
$line_height -$line_height rlineto
-$line_height -$line_height rlineto
EOF
      }
    }
    
    print <<EOF;
stroke
EOF
  }
  
  my $key = $f->key;

  if ($shade_left_keys =~ /\b$key\b/) {
    # clip to the size of the feature so that we can do shading easily
    print <<EOF;
gsave

newpath
$x_start $y_top M
$x_start $y_bottom L
$x_end $y_bottom L
$x_end $y_top L
closepath
clip
EOF

    my $number_of_shade_lines = ($x_end - $x_start)/$shade_spacing;

    for (my $i = -10 ; $i < $number_of_shade_lines ; ++$i) {
      
      my $shade_x_start = $x_start + $shade_spacing *  $i;
      my $shade_x_end = $x_start + $feature_height + $shade_spacing * $i;
      
      print <<EOF;
newpath
$shade_x_start $y_top M
$shade_x_end $y_bottom L
closepath
stroke
EOF

    }

    print <<EOF;
grestore
EOF
  }

  
  if ($shade_right_keys =~ /\b$key\b/) {
    # clip to the size of the feature so that we can do shading easily
    print <<EOF;
gsave

newpath
$x_start $y_top M
$x_start $y_bottom L
$x_end $y_bottom L
$x_end $y_top L
closepath
clip
EOF

    my $number_of_shade_lines = ($x_end - $x_start)/$shade_spacing;

    for (my $i = -10 ; $i < $number_of_shade_lines ; ++$i) {

      my $shade_x_start =
      $x_start + $feature_height + $shade_spacing *  $i;
      my $shade_x_end = $x_start + $shade_spacing * $i;
      
      print <<EOF;
newpath
$shade_x_start $y_top M
$shade_x_end $y_bottom L
closepath
stroke
EOF
      
    }

    print <<EOF;
grestore
EOF
  }


  if (defined $right_exon) {
    # draw connecting lines (an inverted V) between the exons
    my $left_exon = $exon;

    my $next_exon_start = 
     ($right_exon->start - $line_start_base) / 
     $line_base_length * $line_length +
     $border_width;
    
    my $feature_mid_point = ($y_top + $y_bottom) / 2;
    
    my $intron_mid_point = ($next_exon_start + $x_end) / 2;
    
    print <<EOF;
0 0 0 setrgbcolor
newpath
$x_end $feature_mid_point moveto
$intron_mid_point $y_top lineto
$next_exon_start $feature_mid_point lineto
stroke
EOF

  }

}

sub draw_label
{
  my ($complement_flag, $dna_line_feature, $gene_name, 
      $line_start_base, $scale_line_y_pos, @exon_ranges) = @_;

  # get the range of the feature so we can place the label

  my $feature_min = 99999999999;
  my $feature_max = -99999999999;

  for my $exon (@exon_ranges) {
    if ($exon->start < $feature_min) {
      $feature_min = $exon->start;
    }

    if ($exon->end > $feature_max) {
      $feature_max = $exon->end;
    }
    if ($exon->end < $feature_min) {
      $feature_min = $exon->end;
    }

    if ($exon->start > $feature_max) {
      $feature_max = $exon->start;
    }
  }
    
  my $feature_mid =
    (($feature_min + $feature_max) / 2  - $line_start_base) / 
    $line_base_length * $line_length + $border_width;;

  if ($feature_mid >= $border_width && 
      $feature_mid <= $line_length + $border_width) {
    if ($dna_line_feature || $horizontal_labels) {
      my $text_y;
      if ($dna_line_feature) {
        if ($dna_feature_labels eq "on_dna") {
          $text_y = $scale_line_y_pos;
        } else {
          $text_y = $scale_line_y_pos + $feature_height / 2 + 
            $label_feature_distance + $font_size / 2;
        }
      } else {
        if ($complement_flag) {
          $text_y =
            $scale_line_y_pos - $feature_height - $scale_feature_distance -
            $label_feature_distance - $font_size / 2;
        } else {
          $text_y = $scale_line_y_pos + $feature_height + 
            $scale_feature_distance + $label_feature_distance + $font_size / 2;
        }
      }
      print <<EOF;
0 0 0 setrgbcolor
gsave
/Helvetica findfont     $font_size scalefont setfont

$feature_mid $text_y moveto

gsave
newpath
0 0 moveto
($gene_name) false charpath % flattenpath
pathbbox
grestore

neg 2 div exch neg 2 div exch

rmoveto
pop pop
($gene_name)
show
grestore
EOF
    } else {
      if ($complement_flag) {
        my $text_y =
          $scale_line_y_pos - $feature_height - $scale_feature_distance -
          $label_feature_distance;
        print <<EOF;
0 0 0 setrgbcolor
gsave
/Helvetica findfont     $font_size scalefont setfont
gsave
newpath
0 0 moveto
($gene_name) false charpath % flattenpath
pathbbox
grestore

neg exch neg exch

%dup 1 1 roll
%exch dup 1 1 roll exch


$feature_mid $text_y translate

$label_angle rotate
moveto

pop pop

%$feature_mid $text_y moveto

%($gene_name) $label_angle rotate

($gene_name)

show
grestore
EOF
      } else {
        my $text_y = $scale_line_y_pos + $feature_height + 
          $scale_feature_distance + $label_feature_distance;
        print <<EOF;
0 0 0 setrgbcolor
gsave
/Helvetica findfont     $font_size scalefont setfont
$feature_mid $text_y moveto
($gene_name) $label_angle rotate
show
grestore
EOF
      }
    }
  } 
}


sub draw_scale
{
  my ($scale_line_y_pos, $line_start_base, 
      $line_end_base, $line_base_length) = @_;
  
  my $line_x_length = 
    ($line_end_base - $line_start_base) / $line_base_length * $line_length;

  print <<EOF;
0 0 0 setrgbcolor
1.0 setlinewidth
$border_width $scale_line_y_pos M
$line_x_length 0 rlineto
stroke

/Helvetica findfont     $scale_font_size scalefont setfont

EOF

  if ($base_label_position eq "both" ||
      $base_label_position eq "left") {
  print <<EOF;
$border_width $scale_line_y_pos M

gsave
0 0 moveto
($line_start_base) false charpath
pathbbox
grestore

neg 2 div exch neg exch rmoveto

-$scale_number_distance 0 rmoveto

($line_start_base) show
EOF
  }


  if ($base_label_position eq "both" ||
      $base_label_position eq "right") {
  print <<EOF;
$border_width $scale_line_y_pos M
$line_x_length 0 rmoveto
$scale_number_distance 0 rmoveto

gsave
0 0 moveto
($line_end_base) false charpath
pathbbox
grestore

neg 2 div exch pop 0 exch rmoveto

($line_end_base) show

EOF
  }

  for (my $i = $line_start_base - 1 ; $i < $line_end_base ; ++$i) {
    if ($i % $thick_mark_spacing == 0) {
      my $x_pos = 
        ($i - $line_start_base) / $line_base_length * $line_length +
        $border_width;

      my $tick_line_top = $scale_line_y_pos + $thick_tick_line_height / 2;

      print <<EOF;
1.0 setlinewidth
newpath
$x_pos $tick_line_top moveto
0 -$thick_tick_line_height rlineto
stroke
EOF
    }
    if ($i % $tick_mark_spacing == 0) {
      my $x_pos = 
        ($i - $line_start_base) / $line_base_length * $line_length +
        $border_width;

      my $tick_line_top = $scale_line_y_pos + $tick_line_height / 2;

      print <<EOF;
0.3 setlinewidth
newpath
$x_pos $tick_line_top moveto
0 -$tick_line_height rlineto
stroke
EOF
    }
  }
}

my $seq_file = shift @ARGV;

my $seq_stream = AnT::Seqfactory->make (-file => $seq_file);

my $entry = $seq_stream->next_seq ();

if (!defined $entry) {
  die "couldn't read an entry from $seq_file\n";
}


print <<EOF;
%!PS-Adobe
initgraphics
0.50 setlinewidth
/M {moveto} def /L {lineto} def
% newpath
% $border_width $border_width moveto
% $line_length 0 rlineto
% 0 $page_height_less_border rlineto
% -$line_length 0 rlineto
% closepath
% clip
EOF

my $sequence_length = $entry->length ();

# print out each line
for (my $line = 0 ; ; ++$line) {

  my $line_start_base = $line * $line_base_length + 1;

  if ($line_start_base > $sequence_length) {
    last;
  }

  my $line_end_base = $sequence_length;

  if ($line_start_base + $line_base_length < $line_end_base) {
    $line_end_base = $line_start_base + $line_base_length;
  } 

  my $scale_line_y_pos = $page_height - $border_width - $line * $line_spacing -
    $line_spacing / 2;

#  print "line: $line_start_base to $line_end_base\n";

  draw_scale ($scale_line_y_pos, $line_start_base, $line_end_base, 
              $line_base_length);

  my $line_base_range = new AnT::Range (-start => $line_start_base,
                                        -end => $line_end_base);

  my @features = get_features_in_range ($entry, $line_base_range);

# loop over the features on this line
  for my $f (@features) {

    my $key = $f->key;

    if ($protein_features !~ /\b$key\b/ && 
        $dna_line_features !~ /\b$key\b/) {
      # go straight to the next feature
      next;
    }

    if ($filter_in ne "") {
      my ($filter_in_key, $filter_in_qualifier, $filter_in_pattern) =
        split (/:/, $filter_in);

      if ($filter_in_key =~ /\b$key\b/) {

        my @temp_values = $f->qvalues ($filter_in_qualifier);
        
        if (@temp_values) {
          if ($temp_values[0] !~ /$filter_in_pattern/) {
            next;
          }
        }
      }
    }
      
    my $gene_name;

    my $label = $f->label;
      
    if (defined $label) {
      if ($label eq "*") {
        $gene_name = "";
      } else {
        $gene_name = $label;
      }
    } else {
      $gene_name = $f->gene;

      if (!defined $gene_name) {
        $gene_name = $key;
      }
    }

    print <<EOF;
gsave

newpath
$border_width $border_width moveto
$line_length 0 rlineto
0 $page_height_less_border rlineto
-$line_length 0 rlineto
closepath
clip
EOF

    my $dna_line_feature = 0;

    if ($dna_line_features =~ /\b$key\b/) {
      $dna_line_feature = 1;
    }

    $gene_name =~ s/"(.*)"/$1/;

#    print "feature: $$f{key}   $$f{location}\n";
    
    if ($f->strand == 0) {
      warn "ignoring ", $f->key, " at location ", $f->location, "\n";
      next;
    }

    my ($complement_flag, @exon_ranges) = ($f->strand == -1, $f->ranges);

    for (my $exon_index = @exon_ranges -1 ; $exon_index >= 0 ; --$exon_index) {
      draw_exon ($complement_flag, $dna_line_feature,
                 $scale_line_y_pos, $line_start_base,
                 $exon_ranges[$exon_index], 
                 ($exon_index < @exon_ranges - 1) ?
                   $exon_ranges[$exon_index+1] :
                   undef,
                 $exon_index, $f);
    }

    print <<EOF;
grestore
EOF

    if ($dna_line_feature) {
      if ($dna_feature_labels eq "none" || $dna_feature_labels eq "no") {
        next;
      }
    } else {
      if ($protein_feature_labels eq "none" || 
          $protein_feature_labels eq "no") {
        next;
      }
    }

    draw_label ($complement_flag, $dna_line_feature,
                $gene_name, $line_start_base,
                $scale_line_y_pos, @exon_ranges);
  }
}

print <<EOF;
showpage
EOF
