package Bio::Graphics::Browser2::PadAlignment;

use strict;
use Bio::Graphics::Browser2::Markup;
use constant DEBUG=>0;
use Data::Dumper;

=head1 NAME

Bio::Graphics::Browser2::PadAlignment - Insert pads into a multiple alignment

=head1 VERSION (CVS-info)

 $RCSfile: PadAlignment.pm,v $
 $Revision$
 $Author: lstein $
 $Date: 2008-09-09 20:53:36 $

=head1 SYNOPSIS

 use Bio::Graphics::Browser2::PadAlignment;

 my @dnas = (
	     dna1 =>'FFFFgatcGATCgatcGATCgatcGATCgatcGBATCgatcGATCatcGATCgatcGATCgatcGATCgatcGATgatcGATCgatcNNNNGATC',
	     dna2 =>'FFgatcGATCGATCgatcNNGATCgatcGATCgatcGATCgatcGATCgatcGATCtcGATBCgatcGATCatcGATCgatcNNNNGATCFFFF',
	     dna3 =>'FFFFgatcGATCgatcGATCgatcGATCgatcBBBGATCgatcGATCatcGATCBgatcGATCgatcGATCgatcGATgatcGATCgatcNNNNGATCFF',
	     dna4 =>'ZZFFFFgatcGATCgatcGATCgatcGATCgatc',
	     dna5 =>'ATBGGATtcttttttt',
	   );

 #                   target  st  en  tst ten
 my @alignments = ([ 'dna2', 4,  11,  2, 9      ],
 		   [ 'dna2', 16, 23,  10, 17    ],
		   [ 'dna2', 24, 32,  20, 28    ],
		   [ 'dna2', 34, 44,  29, 39    ],
		   [ 'dna2', 45, 59,  41, 55    ],
		   [ 'dna2', 62, 66,  56, 60    ],
		   [ 'dna2', 67, 74,  62, 69    ],
		   [ 'dna2', 76, 86,  71, 81    ],
		   [ 'dna2', 91, 94,  86, 89    ],

		   [ 'dna3',  4,  31, 4, 31    ],
		   [ 'dna3',  33, 33, 34, 34    ],
		   [ 'dna3',  34, 51, 36, 53    ],
		   [ 'dna3',  52, 86, 55, 89    ],
		   [ 'dna3',  91, 94, 94, 97    ],

		   [ 'dna4',  0,  31, 2,  33    ],

		   [ 'dna5',  17, 18, 0, 1    ],
		   [ 'dna5',  41, 43, 4,   6    ],
		   [ 'dna5',  85, 86, 7,   8    ],
		  );


 my $align = Bio::Graphics::Browser2::PadAlignment->new(\@dnas,\@alignments);

 my @padded = $align->padded_sequences;
 print join "\n",@padded,"\n";
 # ..FFFFgatcGATCgatcGATCgatc--GATCgatcG-B-ATCgatcGATC-atcGATC-gatcGATCgatcGAT-CgatcGATgatcGATCgatcNNNNGATC....
 # ....FFgatcGATC----GATCgatcNNGATCgatcG---ATCgatcGATCgatcGATC-gatcGATC--tcGATBCgatcGATCatcGATCgatcNNNNGATCFFFF
 # ..FFFFgatcGATCgatcGATCgatc--GATCgatcBBBGATCgatcGATC-atcGATCBgatcGATCgatcGAT-CgatcGATgatcGATCgatcNNNNGATCFF..
 # ZZFFFFgatcGATCgatcGATCgatc--GATCgatc........................................................................
 # ...................AT-----------------BG-------GAT--------------------------------------------tcttttttt.....

 my $pretty = $align->alignment;
 print $pretty,"\n";
 # dna1   1 ..FFFFgatc GATCgatcGA TCgatc--GA TCgatcG-B- ATCgatcGAT C-atcGATC- gatcGATCga tcGAT-Cgat
 # dna2   1 ....FFgatc GATC----GA TCgatcNNGA TCgatcG--- ATCgatcGAT CgatcGATC- gatcGATC-- tcGATBCgat
 # dna3   1 ..FFFFgatc GATCgatcGA TCgatc--GA TCgatcBBBG ATCgatcGAT C-atcGATCB gatcGATCga tcGAT-Cgat
 # dna4   1 ZZFFFFgatc GATCgatcGA TCgatc--GA TCgatc.... .......... .......... .......... ..........
 # dna5   1 .......... .........A T--------- ---------- -----BGGAT ---------- ---------- ----------

 # dna1  72 cGATgatcGA TCgatcNNNN GATC....
 # dna2  67 cGATCatcGA TCgatcNNNN GATCFFFF
 # dna3  75 cGATgatcGA TCgatcNNNN GATCFF..
 # dna4  35 .......... .......... ........
 # dna5   8 ---------- ----tctttt ttt.....

=head1 DESCRIPTION

This is a utility module for pretty-printing the type of alignment
that comes out of gbrowse, namely a multiple alignment in which each
target is aligned to a reference genome without explicit pads or
other spaces.

For speed and ease of use, the module does not use form Bio::SeqI
objects, but raw strings and alignment data structures.  This may
change.

This module does B<not> perform multiple alignments!  It merely
pretty-prints them!

=head2 METHODS

This section describes the methods used by this class.

=over 4

=item $aligner = Bio::Graphics::Browser2::PadAlignment->new(\@sequences,\@alignments)

Create a new aligner.  The two arguments are \@sequences, an array ref
to the list of sequences to be aligned, and \@alignments, an array ref
describing how the sequences are to be aligned.

\@sequences should have the following structure:

  [ name1 => $sequence1,
    name2 => $sequence2,
    name3 => $sequence3 ]

The sequences will be displayed in top to bottom order in the order
provided.  The first sequence in the list is special because it is the
reference sequence.  All alignments are relative to it.

\@alignments should have the following structure:

 [ [ target1, $start1, $end1, $tstart1, $tend1 ],
   [ target1, $start2, $end2, $tstart2, $tend2 ],
   ...
  ]

Each element of @alignments is an arrayref with five elements.  The
first element is the name of the target sequence, which must be one of
the named sequences given in @sequences.  The second and third
elements are the start and stop position of the aligned target segment
relative to the reference sequence, "name1" in the example given
above.  The fourth and fifth elements are the start and stop position
of the aligned target segment in the coordinate space of the target.

Example:

  @dnas = ('dna1' => 'ccccccaaaaaatttt',
	   'dna2' => 'aaaaaa');
  @alignment = ( ['dna2', 6, 11, 0, 5 ]);

Positions 0 to 5 of "dna2" align to positions 6-11 of "dna1".

Note that sequence positions are zero based.  This may change.

=item @lines = $aligner->padded_sequences

This inserts pads into the sequences and returns them as a list of
strings in the order specified in new().  In a scalar context, this
method will return a hashref in which the keys are the sequence names
and the values are their padded strings.

=item $map = $aligner->gap_map

This returns an arrayref indicating the position of each base in the
gapped reference sequence.  The indexes are base positions, and the
element values are their positions in the reference sequence as
returned by padded_sequences().

Note that the gap map only provides coordinate mapping for the
reference sequence.  For an alternative implementation that provides
gap maps for each of the targets (at the cost of speed and memory
efficiency) see the section after __END__ in the source file for this
module.

=item $align_string = $aligner->alignment(\%origins [,\%options])

This method returns a pretty-printed string of the aligned sequences.
You may provide a hashref of sequence origins in order to control the
numbers printed next to each line of the alignment.  The keys of the
%origins hashref are the names of the sequences, and the values are
the coordinate to be assigned to the first base of the sequence.  Use
a negative number if you wish to indicate that the sequence has been
reverse complemented (the negative number should indicate the
coordinate of the first base in the provided sequence).

An optional second argument, if present, contains a hash reference to
a set of option=>value pairs.  Three options are recognized:

   show_mismatches      0|1      if true, highlight mismatches in pink
   show_matches		0|1	 if true, hightligt matches in yellow
   color_code_proteins  0|1      if true, highlight amino acids thus:
                                        Acidic amino acids in red
                                        Basic amino acids in blue
                                        Hydrophobic amino acids in grey
                                        Polar amino acids in yellow
   flip                 0|1      if true, reverse complement the whole alignment

=back


=head1 SEE ALSO

L<Bio::Graphics::Panel>,
L<Bio::Graphics::Glyph>,
L<Bio::Graphics::Feature>,
L<Bio::Graphics::FeatureFile>,
L<Bio::Graphics::Browser>,
L<Bio::Graphics::Browser2::Plugin>

=head1 AUTHOR

Lincoln Stein E<lt>lstein@cshl.orgE<gt>.

Copyright (c) 2002 Cold Spring Harbor Laboratory

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself.  See DISCLAIMER.txt for
disclaimers of warranty.

=cut

# A specific package for padding multiple alignments into an ASCII text string.
# It is designed for cases in which all alignments are against
# a single reference sequence, such as ESTs against a genome.  The reference
# sequence must be the first one provided.

# for efficiency, we use zero-based coordinates throughout.

# IMPORTANT NOTE: see the section after __END__ for a slightly
# different implementation which keeps a separate gap map for each 
# sequence in the alignment

# define the types of amino acids -- this was done by an undergrad and is subject to change
# modified according to: http://www.ann.com.au/MedSci/amino.htm, method 1

my %aa_type = (
	       K=> "basic_aa",
	       R=> "basic_aa",
	       H=> "basic_aa",
	       S=> "polar_aa",
	       T=> "polar_aa",
	       N=> "polar_aa",
	       Q=> "polar_aa",
	       C=> "polar_aa",
	       Y=> "polar_aa",
	       D=> "acidic_aa",
	       E=> "acidic_aa",
	       G=> "npolar_aa",
	       A=> "npolar_aa",
	       V=> "npolar_aa",
	       L=> "npolar_aa",
	       I=> "npolar_aa",
	       P=> "npolar_aa",
	       M=> "npolar_aa",
	       F=> "npolar_aa",
	       C=> "npolar_aa",
	       W=> "npolar_aa",
	       X=> "special_aa",
	       "*" => "special_aa"
	      );

sub new {
  my $class = shift;
  my $dnas   = shift;  # array ref of DNAs in the order in which they will be printed
                       # in format [ [name1=>dna1],[name2=>dna2]...]
  my $aligns = shift;  # array ref of alignments in format [ [targetname,srcstart,srcend,targetstart,targetend] ]

  # remap data structures
  my $count = 0;
  my (@dnas,%dnas);
  while (my($name,$dna) = splice(@$dnas,0,2)) {
    $dnas{$name} = $count++;
    push @dnas,$dna;
  }
  return bless {
		names  => \%dnas,
		dnas   => \@dnas,
		aligns => $aligns
		};
}

# return a hashref in which the keys are the
# dna names and the values are the padded strings
sub padded_sequences {
  my $self = shift;
  my @lines;

  my @dnas = @{$self->{dnas}};

  # initialize top line with src sequences
  $lines[0] = $dnas[0];
  my $len = length($lines[0]); 
  for (my $i = 1; $i < @dnas; $i++) {
    $lines[$i]  = '-' x $len;
  }

  # running total of number of gaps, indexed by position on source
  my @gap_map = 0..length($dnas[0])-1;

  # place where DNA[$i] left off
  my @last_end;
  my @added = (length($dnas[0])-1);

  # leader sequence to add back after gapping
  my @leader = (0);

  # alignments must be sorted according to target
  foreach (sort {$a->[0] cmp $b->[0]
		   || $a->[1] <=> $b->[1]
		 } @{$self->{aligns}}) {
    my ($targ, $start, $end,
	$tstart,$tend) = @{$_};

    defined ($targ = $self->{names}{$targ}) or next;

    warn "$start $end $tstart $tend\n" if DEBUG;

    my $src = 0;   # first DNA is the reference

    # position in src coordinates where we last stopped
    my $last_src    = $last_end[$targ][$src]  ||= -1;

    # position in target coordinates where we last stopped
    my $last_target = $last_end[$targ][$targ] ||= -1;

    if ($last_src >= 0 && $last_target >= 0) {
      # This section adds the unaligned region between the last place that
      # we stopped and the current alignment
      my $gap       = 0;
      my $tpos      = $last_target+ 1;
      my $spos      = $gap_map[$last_src]+1;
      warn "last_src=$last_src, spos=$spos" if DEBUG;
      my $deficit   = ($tstart-$last_target) - ($gap_map[$start] - $gap_map[$last_src]);
      warn "add $deficit gaps" if DEBUG;
      if ($deficit > 0) {
	for (my $i=0; $i<@lines; $i++) {
	  eval {substr($lines[$i],$spos,0) = '-'x$deficit};
	}
	@gap_map[$start..$#gap_map] = map {$_+$deficit} @gap_map[$start..$#gap_map];
      }
      while ($tpos < $tstart) {
	eval { substr($lines[$targ],$spos++,1) = substr($dnas[$targ],$tpos++,1) };
      }
    }

    else {  # remember to add the extra stuff at beginning
      $leader[$src]  = $start  if !defined $leader[$src]  || $start < $leader[$src];
      $leader[$targ] = $tstart if !defined $leader[$targ] || $tstart< $leader[$targ];
    }

    # insert the aligned bit now
    for (my $pos = $start; $pos <= $end; $pos++) {
      my $gap_pos = $gap_map[$pos];
      defined $gap_pos or next;
      warn "inserting $gap_pos with ",substr($dnas[$targ],$tstart,1),"\n" if DEBUG;
      eval {substr($lines[$targ],$gap_pos,1) = substr($dnas[$targ],$tstart++,1) };
    }

    $last_end[$targ][$src]  = $end;
    $last_end[$targ][$targ] = $tend;

    $added[$targ] = $tend;

    warn join("\n",@lines),"\n\n" if DEBUG;
  }

  # take care of the extra stuff at the end
  for (my $i=1; $i < @dnas; $i++) {
    my $last_bit = length($dnas[$i]) - $added[$i];
    next unless defined $gap_map[$last_end[$i][0]];
    local $^W = 0; # prevent uninit variable warnings
    eval {substr($lines[$i],$gap_map[$last_end[$i][0]]+1,$last_bit)
	    = substr($dnas[$i],$added[$i]+1,$last_bit) };
  }

  # take care of the extra unaligned stuff at the beginning
  my $max = 0;
  for (my $i=0; $i < @dnas; $i++) {
    $leader[$i] ||= 0;      # to prevent uninit variable warnings
    next unless $leader[$i];
    my ($leading_gaps) = $lines[$i] =~ /^(-+)/;
    my $leading_pads   = length($leading_gaps||'');

    warn "\$leader[$i] = $leader[$i], \$leading_pads = $leading_pads\n" if DEBUG;

    my $insert_length = $leading_pads >= $leader[$i] ? $leader[$i] : $leading_pads;
    my $append_length = $leading_pads >= $leader[$i] ? 0           : $leader[$i]-$leading_pads;
    warn "insert length = $insert_length, append_length=$append_length\n" if DEBUG;

    if ($insert_length > 0) {
      substr($lines[$i],$leading_pads-$insert_length,$insert_length) =
	substr($dnas[$i],$leader[$i]-$insert_length,$insert_length);
      $leader[$i] -= $insert_length;
    }
    if ($append_length > 0) {
      substr($lines[$i],0,0) = $i>0 ? substr($dnas[$i],0,$append_length) : '-'x$append_length;
    }

    $max = $append_length if $append_length > $max;
  }
  warn "\n" if DEBUG;
  warn join("\n",@lines),"\n\n" if DEBUG;

  warn "finished adding stuff for everything but reference sequence\n" if DEBUG;

  for (my $i=0; $i<@dnas; $i++) {
    warn "i = $i, max = $max, leader = $leader[$i]\n" if DEBUG;
    my $delta = $max - $leader[$i];
    next unless $delta > 0;
    substr($lines[$i],0,0) = '-'x$delta;
  }

  warn join("\n",@lines),"\n\n" if DEBUG;

  # change starts and ends to . characters
  $max = 0;
  foreach (@lines) {
    $max = length if $max < length;
  }
  foreach (@lines) {
    my $deficit = $max - length;
    s/^(-+)/'.'x length $1/e;
    s/(-+)$/'.'x length $1/e;
    $_ .= '.' x $deficit if $deficit > 0;
  }

  $self->{gaps} = \@gap_map;

  return @lines if wantarray;

  my %names = reverse %{$self->{names}};
  my %result = map {$names{$_} => $lines[$_]} 0..$#lines;
  return \%result;
}

sub gap_map {
  my $self  = shift;
  $self->padded_sequences unless $self->{gaps};
  my $map   = $self->{gaps};
  return $map;
}

sub alignment {
  my $self            = shift;
  my $origins         = shift;
  my $options         = shift || {};

  my $show_mismatches = $options->{show_mismatches};
  my $show_matches = $options->{show_matches};
  my $show_similarities = $options->{show_similarities};
  my $color_code_proteins = $options->{color_code_proteins};
  warn "color code = $color_code_proteins" if DEBUG;
  my $flip            = $options->{flip};

  my @lines = $self->padded_sequences;
  my %names = reverse %{$self->{names}};  # index to name

  $origins ||= {};

  foreach (values %names) {
    $origins->{$_} = 1 unless defined $origins->{$_};
  }

  my $longest_name = 0;
  foreach (values %names) {
    my $offset    = $origins->{$_};
    my $length    = length($_);
    $length      += 2 if $offset < 0;
    $longest_name = $length if $length > $longest_name;
  }
  my $longest_line = 0;
  for (my $i=0; $i<@lines; $i++) {
    my $offset    = abs($origins->{$names{$i}});
    $longest_line = length($self->{dnas}[$i])+$offset 
       if (length($self->{dnas}[$i])+$offset > $longest_line);
  }

  $longest_line = length $longest_line;  # looks like an error but isn't

  # if flip is set, then we do amazing things to reorganize the display!
  if ($flip) {
    for (my $i = 0; $i < @lines; $i++) {
      $lines[$i] = reverse $lines[$i];
      $lines[$i] =~ tr/gatcGATC/ctagCTAG/;
      my $name   = $names{$i};
      $origins->{$name} *= -1;
    }
  }
  
  # use markup to insert word and line breaks
  my $markup = Bio::Graphics::Browser2::Markup->new;
  $markup->add_style(space    => ' ');
  $markup->add_style(newline  => "\n");
  $markup->add_style(mismatch => "BGCOLOR pink");
  $markup->add_style(match => "BGCOLOR darkorange");
  $markup->add_style(conserved => "BGCOLOR tan");

  # Styles for printing protein alignments
  $markup->add_style(acidic_aa => "BGCOLOR lightgreen");
  $markup->add_style(basic_aa => "BGCOLOR lightskyblue");
  $markup->add_style(npolar_aa => "BGCOLOR lightgrey");
  $markup->add_style(polar_aa => "BGCOLOR burlywood");
  $markup->add_style(special_aa => "BGCOLOR red");

  # add word and line breaks
  for (my $i=0; $i < @lines; $i++) {
    my $pad = \$lines[$i];
    my @markup;
    for (my $j=0; $j < length $$pad; $j += 10) {
      push (@markup,[$j % 80 ? 'space':'newline',
                     $j => $j]);
    }
    $markup->markup($pad,\@markup);
  }

  my (@padded, @labels, @fixed);

  for (my $i = 0; $i < @lines; $i++) {
    my @segments = split "\n",$lines[$i];
    for (my $j = 0; $j < @segments; $j++) {
      $padded[$j][$i] = $segments[$j];
      $fixed[$j][$i] = $segments[$j];
    }
    my $origin  = $origins->{$names{$i}};
    $labels[$i] = $origin                                      if $origin >= 0;
    $labels[$i] = length($self->{dnas}[$i]) + abs($origin) - 1 if $origin <  0;
  }

  my $result;
  my @length;

  #$i: number of blocks; $j: number of sequences
  for (my $i = 0; $i < @padded; $i++) {
    for (my $j = 0; $j < @{$padded[$i]}; $j++) {

      next unless $padded[$i][$j];
      my $origin = $origins->{$names{$j}};
      my $offset = $padded[$i][$j] =~ tr/. -/. -/;
      my $skipit = $offset == length($padded[$i][$j]); 
      my @markup;

      #warn "Block ", $i, "\tsequence ", $j, "\t", $origin, "\t", $offset, "\t", $skipit, "\n";

      if ($color_code_proteins){
        if ($j==0) {                            # colouring reference seq
          for(my $q=0; $q<length $padded[$i][$j]; $q++) {
            my $refPos = substr($padded[$i][$j],$q,1);
            next if $refPos =~ /^[.\s-]$/;                               # move on if not amino acid
            push(@markup,[$aa_type{$refPos},$q=>$q+1]);
          } # end for
        }

	else {
          my @markup;
          for (my $r=0; $r<length $padded[$i][$j]; $r++) {
            my $targ = substr($padded[$i][$j],$r,1);
            next if $targ =~  /^[.\s-]$/;
            push(@markup,[$aa_type{$targ}, $r => $r+1]);
          }
        }
      }

      elsif ($show_mismatches) {
        if ($j>0) {
          for (my $r=0; $r<length $padded[$i][$j]; $r++) {
            my $targ = substr($padded[$i][$j],$r,1);
            next if $targ =~  /^[.\s-]$/;
            my $source = substr($padded[$i][0],$r,1);
            next if $source=~ /^[.\s-]$/;

            push(@markup,['mismatch',$r => $r+1])
              if (lc($source) ne lc($targ)); 

          }
        }
      }

      elsif ($show_matches) {
        for (my $r=0; $r<length $padded[$i][$j]; $r++) {

          my $targ = substr($padded[$i][$j],$r,1);
          my $identical = 1;
          my $conserved = 1;

          for (my $m=0; $m<@{$fixed[$i]}; $m++){
            my $source = substr($fixed[$i][$m],$r,1);

            if(($source =~ /[.\s-]/)||($targ =~ /[.\s-]/)){
              $identical = undef;
            } elsif (lc($source) ne lc($targ)){
              $identical = undef;
            }

            if(($source =~ /[.\s-]/)||($targ =~ /[.\s-]/)){
              $conserved = undef;
            }elsif (($aa_type{$source}) ne ($aa_type{$targ})){
              $conserved = undef;
            }
          }
          push(@markup,['match',$r => $r+1]) if ($identical);
          push(@markup,['conserved',$r => $r+1]) if ($conserved);
        }
      }

      elsif ($show_similarities) { #highligt resides same to the reference protein
        if ($j == 0){
          for (my $r=0; $r<length $padded[$i][$j]; $r++) {
            my $identical = undef;
            my $conserved = undef;

            my $targ = substr($padded[$i][$j],$r,1);
            next if ($targ =~ /^[.\-]$/);

            for (my $m=1; $m<@{$fixed[$i]}; $m++){
              my $source = substr($fixed[$i][$m],$r,1);
              next if ($source =~ /^[.\-]$/);

              if (($source !~ /^[.\s-]$/) && ($targ !~ /^[.\s-]$/) && (lc($source) eq lc($targ))){
                $identical = 1;
              }

              if (($source !~ /^[.\s-]$/) && ($targ !~ /^[.\s-]$/) && 
                  (lc($aa_type{$source}) eq lc($aa_type{$targ})) &&
                  (lc($source) ne lc($targ))
                 ){
                $conserved = 1;
              }
            }
            push(@markup,['conserved',$r => $r+1]) if ($conserved);
            push(@markup,['match',$r => $r+1]) if ($identical);
          }
        }

	else {
          for (my $r=0; $r<length $padded[$i][$j]; $r++) {

            my $targ = substr($padded[$i][$j],$r,1);
            my $identical = undef;
            my $conserved = undef;

            my $source = substr($fixed[$i][0],$r,1);
            next if ($source =~ /^[.\-]$/);

            if (($source !~ /[.\s-]/) && (lc($source) eq lc($targ))){
              $identical = 1;
            }

            if (($source !~ /^[.\s-]$/) && ($targ !~ /^[.\s-]$/) && 
                (lc($aa_type{$source}) eq lc($aa_type{$targ})) &&
                (lc($source) ne lc($targ))
               ){
              $conserved = 1;
            }

            push(@markup,['conserved',$r => $r+1]) if ($conserved);
            push(@markup,['match',$r => $r+1]) if ($identical);
          }
        }
      }

      $length[$i][$j] = length $padded[$i][$j];
      $markup->markup(\$padded[$i][$j],\@markup) if @markup;

      my $l = $longest_name;
      $result .= $skipit ? ""
                       : sprintf ("\%${l}s \%${longest_line}d %s\n",
                                  $origin < 0 ? "($names{$j})"
                                                : $names{$j},
                                    $labels[$j],$padded[$i][$j]);

      $labels[$j] += $length[$i][$j] - $offset  if $origin >= 0;
      $labels[$j] -= $length[$i][$j] - $offset  if $origin < 0;
    }
    $result .= "\n"; # unless $result && $result =~ /^[.\s]+$/;  # skip completely empty lines
  }	
  return $result;
}

1;

__END__

use constant SRC   => 0;
use constant TARG  => 1;

use constant REF   => 0;
use constant START => 1;
use constant END   => 2;

my @dnas = (
#               10        20          30          40         50         60         70         80        90
#      012345678901234567890123  456789012 3 45678901234 5678901 234567890123456 78901234 56789012345678901234
      'FFFFgatcGATCgatcGATCgatc--GATCgatcG-B-ATCgatcGATC-atcGATC-gatcGATCgatcGAT-CgatcGAT-gatcGATCgatcNNNNGATC',
        'FFgatcGATC----GATCgatcNNGATCgatcG---ATCgatcGATCgatcGATC-gatcGATC--tcGATBCgatcGATC-atcGATCgatcNNNNGATCFFFF',
#        0123456789    0123456789012345678   9012345678901234567 89012345  678901234567890 12345678901234567890123
      'FFFFgatcGATCgatcGATCgatc--GATCgatcBBBGATCgatcGATC-atcGATCBgatcGATCgatcGAT-CgatcGAT-gatcGATCgatcNNNNGATCFF',
#      012345678901234567890123  45678901234567890123456 78901234567890123456789 01234567 8901234567890123456789
#               10        20          30        40         50        60         70         80        90
    'ZZFFFFgatcGATCgatcGATCgatc--GATCgatc',
#    01234567890123456789012345  67890123
#             10        20          30
	               'AT-----------------BG-------GAT---------------------------------------------tcttttttt',
#                       01                 23       456                                             78
	   );

#                  ref st  en    tar tst ten
my @alignments = ([ [0, 4,  11],  [1,  2, 9]      ],
		  [ [0, 16, 23],  [1,  10, 17]    ],
		  [ [0, 24, 32],  [1,  20, 28]    ],
		  [ [0, 34, 44],  [1,  29, 39]    ],
		  [ [0, 45, 59],  [1,  41, 55]    ],
		  [ [0, 62, 66],  [1,  56, 60]    ],
		  [ [0, 67, 74],  [1,  62, 69]    ],
		  [ [0, 76, 86],  [1,  71, 81]    ],
		  [ [0, 91, 94],  [1,  86, 89]    ],

		  [ [0,  4, 31],  [2,   4, 31]    ],
		  [ [0,  33, 33], [2,  34, 34]    ],
		  [ [0,  34, 51], [2,  36, 53]    ],
		  [ [0,  52, 86], [2,  55, 89]    ],
		  [ [0,  91, 94], [2,  94, 97]    ],

		  [ [0,  0,  31], [3,  2,  33]    ],

		  [ [0,  17, 18], [4,  00, 01]    ],
		  [ [2,  34, 35], [4,  2,   3]    ],
		  [ [0,  41, 43], [4,  4,   6]    ],
		  [ [0,  85, 86], [4,  7,   8]    ],

		 );

foreach (@dnas) { s/-//g };

my @lines;

# initialize top line with src sequences
$lines[0] = $dnas[0];
my $len = length($lines[0]); 
for (my $i = 1; $i < @dnas; $i++) {
  $lines[$i]  = '-' x $len;
}

# running total of number of gaps, indexed by position on source
my @gap_map;
foreach (@dnas) { 
  push @gap_map,[0..length($_)-1]
}

# place where DNA[$i] left off
my @last_end;
my @added = (length $dnas[0]-1);

# alignments must be sorted according to target
foreach (sort {$a->[TARG][REF] <=> $b->[TARG][REF]
		 || $a->[TARG][START] <=> $b->[TARG][START]
	       } @alignments) {
  my ($src, $start, $end)  = @{$_->[SRC]};
  my ($targ,$tstart,$tend) = @{$_->[TARG]};

  # position in src coordinates where we last stopped
  my $last_src    = $last_end[$targ][$src]  ||= -1;

  # position in target coordinates where we last stopped
  my $last_target = $last_end[$targ][$targ] ||= -1;

  # This section adds the unaligned region between the last place that
  # we stopped and the current alignment
  my $gap       = 0;
  for (my $targ_pos=$tstart-1, my $j=$start-1; $targ_pos > $last_target; $targ_pos--, $j--) {
    if ($j > $last_src) { # still room
      my $pos = $gap_map[$src][$j];
      substr($lines[$targ],$pos,1) = substr($dnas[$targ],$targ_pos,1);
    }
    else {  # we've overrun -- start gapping above
      my $pos = $gap_map[$src][$start];
      for (my $i=0; $i<@lines; $i++) {
	substr($lines[$i],$pos,0) = '-' unless $i==$targ;  # gap all segments
      }
      substr($lines[$targ],$pos+$gap++,0) = substr($dnas[$targ],$targ_pos,1);
    }
  }
  if ($gap > 0) {
    for (my $i=0; $i<@lines; $i++) {
      next if $i == $targ;
      for (@{$gap_map[$i]}[$start..$#{$gap_map[$i]}]) { $_ += $gap }  # update gap map
    }
  }

  # insert the aligned bit now
  for (my $pos = $start; $pos <= $end; $pos++) {
    my $gap_pos = $gap_map[$src][$pos];
    substr($lines[$targ],$gap_pos,1) = substr($dnas[$targ],$tstart++,1);
  }

  $last_end[$targ][$src]  = $end;
  $last_end[$targ][$targ] = $tend;

  $added[$targ] = $tend;

  #print join("\n",@lines),"\n\n";
}

# take care of the extra stuff at the end
for (my $i=1; $i < @dnas; $i++) {
  my $last_bit = length($dnas[$i]) - $added[$i];
  substr($lines[$i],$gap_map[0][$last_end[$i][0]]+1,$last_bit)
    = substr($dnas[$i],$added[$i]+1,$last_bit);
}

# change starts and ends to . characters
my $max = 0;
foreach (@lines) {
  $max = length if $max < length;
}
foreach (@lines) {
  my $deficit = $max - length;
  s/^(-+)/'.'x length $1/e;
  s/(-+)$/'.'x length $1/e;
  $_ .= '.' x $deficit if $deficit > 0;
}


print join("\n",@lines),"\n";

# use markup to insert word and line breaks
my $markup = Bio::Graphics::Browser2::Markup->new;
$markup->add_style(space   => ' ');
$markup->add_style(newline => "\n");
for (my $i=0; $i < @lines; $i++) {
  my $pad = \$lines[$i];
  my @markup;
  # add word and line breaks
  for (my $j=0; $j < length $$pad; $j += 10) {
    push (@markup,[$j % 80 ? 'space':'newline',
		   $j => $j]);
  }
  $markup->markup($pad,\@markup);
}

my @padded;
for (my $i = 0; $i < @lines; $i++) {
  my @segments = split "\n",$lines[$i];
  for (my $j = 0; $j < @segments; $j++) {
    $padded[$j][$i] = $segments[$j];
  }
}

my @labels = (1) x @lines;

for (my $i = 0; $i < @padded; $i++) {
  for (my $j = 0; $j < @{$padded[$i]}; $j++) {
    next unless $padded[$i][$j];
    printf ("%5d %s\n",$labels[$j],$padded[$i][$j]);
    my $offset = $padded[$i][$j] =~ tr/. -/. -/;
    $labels[$j] = length($padded[$i][$j]) - $offset + 1;
  }
  print "\n";
}


sub segment {
  my ($str,$start,$stop) = @_;
  my $length = $stop-$start+1;
  return unless $length > 0;
  return substr($str,$start,$length);
}