#!/usr/local/bin/perl

# $Header: /nfs/disk222/yeastpub/Repository/zoo/general/embl_join.pl,v 1.7 1999/11/09 14:16:46 kmr Exp $

# purpose: This script will merge two or more embl entries.  The sequence
# purpose: and features will be merged, the header from the first entry
# purpose: will be used on the output.

use strict;

use Embl;
use EMBL::Parser;

sub usage
{
  die <<EOF;
usage: $0 -h (for help)
   or: $0 [-<num>] [-n] [-o output_file] embl_file1 embl_file2 ...

This script will join two or more EMBL entries in the order given on the
command line.

Options:
-<num>   allows the expected overlap to be given. If the actual overlap is
         different to the number given, a message will be displayed and
         processing will stop.

-n       will force the to entries to be joined without an overlap.

-o file  will send the output to the given file rather than standard output.


example:
    $0 -310 c1919_submission c790_submission > merged_entry.embl

This will produce a new embl file (in merged_entry.embl).  The new file
will contain the header from c1919_submission and the merged sequence and
features from both entries.  If the overlap isn't 310 processing will stop.
EOF
}

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

my $expected_overlap = undef;

if ($ARGV[0] =~ /^-(\d+)$/) {
  $expected_overlap = $1;
  shift;
}

my $no_overlap = 0;

if ($ARGV[0] eq "-n") {
  $no_overlap = 1;
  shift;
}

if ($ARGV[0] eq "-o") {
  if (@ARGV > 1) {
    close STDOUT;
    open STDOUT, ">$ARGV[1]" or die "can not open $ARGV[1]\n";
    shift;
    shift;
  } else {
    usage ();
  }
}

if (@ARGV < 2) {
  usage ();
}



# returns the new sequence and the overlap (in bases)
sub get_merged_sequence
{
  my ($seq_1, $seq_2, $expected_overlap) = ($_[0]->seq, $_[1]->seq, $_[2]);

  my $seq_1_start_pos = 0;

  if (defined $expected_overlap) {
    $seq_1_start_pos = (length $seq_1) - $expected_overlap;
  }

  my $test_base_count = 10;

  for (my $i = $seq_1_start_pos ; $i < length $seq_1 ; ++$i) {

    if ((length $seq_1) - $i > $test_base_count and
        (length $seq_2) > $test_base_count and
        (substr $seq_1, $i, $test_base_count) ne
        (substr $seq_2, 0, $test_base_count)) {
      # failed the quick test
      next;
    }

    # now do a rigorous check
    my $seq_2_substr = substr ($seq_2, 0, length ($seq_1) - $i);
    if (index ($seq_1, $seq_2_substr, $i) == $i) {
      return (substr ($seq_1, 0, $i) . $seq_2,
              (length $seq_1) - $i);
    }
  }

  ($seq_1 . $seq_2, 0);
}


# returns an array containing a complement flag and then a Range for each exon
# ie. for the location "join(complement(1..100),complement(200..400))" the
# return array will be (1, new Range (1,100), new Range (200,400))
sub get_exon_ranges
{
  my $location = shift;

  my @return_vector = ();

  my $complement_flag = 0;

  if ($location =~ s/complement\(([^\)]+)\)/$1/g) {
    $complement_flag = 1;
  }

  $location =~ s/join\(([^\)]*)\)/$1/g;

  my @join_bits = split (',',$location);

  # sort the join ranges into ascending order
  @join_bits =
    sort {
      my ($astart) = $a =~ /^\s*(\d+)/;
      my ($bstart) = $b =~ /^\s*(\d+)/;
      $astart <=> $bstart;
    } @join_bits;

  for (my $i = 0 ; $i < @join_bits ; ++$i) {
    my $exon_range = $join_bits[$i];

    my $start_of_exon;
    my $end_of_exon;

    if ($exon_range =~
        m/<?(\d+)
          (?:\.\.
           >?
           (\d+))?/x) {

      $start_of_exon = $1;
      if (defined $2) {
        $end_of_exon = $2;
      } else {
        $end_of_exon = $1;
      }

      if ($complement_flag) {
        push @return_vector,
          new Range ($end_of_exon, $start_of_exon);
      } else {
        push @return_vector,
          new Range ($start_of_exon, $end_of_exon);
      }
    }
  }

  if ($complement_flag) {
    return ($complement_flag,reverse @return_vector);
  } else {
    return ($complement_flag,@return_vector);
  }
}


sub merge_locations
{
  my ($first_location, $second_location) = @_;

  my ($first_complement_flag, @first_ranges) =
    get_exon_ranges ($first_location);
  my ($second_complement_flag, @second_ranges) =
    get_exon_ranges ($second_location);

  if ($first_complement_flag != $second_complement_flag) {
    # throw a wobbly

    warn "$0: locations are on different strands (does a sequence need ". 
         "reversing?)\n";

    return undef;
  }

  my @merged_ranges = (@first_ranges, @second_ranges);

  # sort the ranges using the start coordinate then the end coordinate
  # of each range
  @merged_ranges = sort {
    my $comp = $a->start <=> $b->end; 
    if ($comp == 0) {
      $a->end <=> $b->end;
    } else {
      $comp;
    }
  } @merged_ranges;

  
  
  # merge overlapping ranges
  for (my $i = @merged_ranges - 1 ; $i >= 1 ; --$i) {
    my $this_range = $merged_ranges[$i];
    my $prev_range = $merged_ranges[$i - 1];

    if ($prev_range->overlapping ($this_range)) {
      my $new_range = new Range ($prev_range->start, $this_range->end);
     
      splice @merged_ranges, $i - 1, 2, $new_range;
    }
  }

  my $return_string = "";

  if ($first_complement_flag) {
    $return_string = "complement(";
  }

  if (@merged_ranges > 1) {
    $return_string .= "join(";
  }

  $return_string .= join ",", map {$_->start . ".." . $_->end} @merged_ranges;

  if ($first_complement_flag) {
    $return_string .= ")";
  }

  if (@merged_ranges > 1) {
    $return_string .= ")";
  }

  return $return_string;
}

sub add_qualifiers
{
  my ($source_feature, $dest_feature) = @_;

  for my $source_qualifier ($source_feature->qualifiers) {
    if (!grep {$_->name eq $source_qualifier->name &&
               lc $_->value eq lc $source_qualifier->value}
        $dest_feature->qualifiers) {
      $dest_feature->addQualifiers ($source_qualifier);
    }
  }
}

sub merge_entries
{
  my ($entry_1, $entry_2) = @_;


  my $sequence_1 = $entry_1->SQ ();
  
  my $sequence_1_length = $sequence_1->length ();
  
  
  my $sequence_2 = $entry_2->SQ ();
  
  my $sequence_2_length = $sequence_2->length ();
  
  
  my ($new_sequence, $overlap);
  
  if ($no_overlap) {
    ($new_sequence, $overlap) = ($sequence_1->seq . $sequence_2->seq, 0);
  
    print STDERR "joined sequences without an overlap\n";
  } else {
    ($new_sequence, $overlap) =
      get_merged_sequence ($sequence_1, $sequence_2, $expected_overlap);
  
    if (defined $expected_overlap and $expected_overlap != $overlap) {
      die "expected overlap ($expected_overlap) is different to the actual " .
          "overlap ($overlap)\n";
    }
  
    if ($overlap == 0) {
      print STDERR "found no overlap\n";
    } else {
      print STDERR "found overlap of $overlap bases\n";
    }
  }
  
  my $overlap_start = $sequence_1->length () - $overlap + 1;
  my $overlap_end = $sequence_1->length ();
  
  if ($overlap != 0) {
    print OVERLAP_TAB <<EOF;
FT   misc_feature    $overlap_start..$overlap_end
FT                   /colour=7
EOF
  }
  
  
  my $overlap_range = new Range ($overlap_start, $overlap_end);
  
  my $next_entry_overlap_range = new Range (1, $overlap);
  
  
  # change the sequence of the first entry
  $sequence_1->seq ($new_sequence);
  
  
  my $ft_1 = $entry_1->FT ();
  
  my $ft_2 = $entry_2->FT ();
  
  my @features_1 = $ft_1->features;
  
  my @features_2 = $ft_2->features;
  
  
  # fix the locations of the features in the second entry by adding
  # $sequence_1_length - $overlap to each base position
  
  my $bases_to_move = $sequence_1_length - $overlap;
  
  for my $feature (@features_2) {
    my $location = $feature->location;
  
    $location =~ s/(\d+)/$1 + $bases_to_move/eg;
  
    $feature->location ($location);
  }
  
  # if the last CDSs of the first have a /partial qualifier then
  # remove them from the entry and remember them so that we can try to
  # join them to the first CDSs of the next entry
  
  my @this_entry_overlap_features = ();
  
  for (my $i = (@features_1 - 1) ; $i >= 0 ; --$i) {
    my $feature = $features_1[$i];
  
    if ($feature->key eq "CDS") {
      if (grep {$_->name =~ /^partial\s*$/} $feature->qualifiers or
          $feature->range ()->overlapping ($overlap_range)) {
        push @this_entry_overlap_features, $feature;
        $ft_1->removeFeature ($i);
      }
    }
  }
  
  # get the first CDS of the second entry.  if they have a /partial
  # qualifier then remove them from the entry and remember them so
  # that we can try to join them to $last_partial_cds
  
  my @next_entry_overlap_features = ();
  
  for (my $i = (@features_2 - 1) ; $i >= 0 ; --$i) {
    my $feature = $features_2[$i];
  
    if ($feature->key eq "CDS") {
      if (grep {$_->name =~ /^partial\s*$/} $feature->qualifiers or
          $feature->range ()->overlapping ($next_entry_overlap_range)) {
        push @next_entry_overlap_features, $feature;
        $ft_2->removeFeature ($i);
      }
    }
  }
  
  
  # try to merge each feature in @this_entry_overlap_features with a
  # feature in @next_entry_overlap_features
  
  this_entry_loop:
    for (my $this_count = (@this_entry_overlap_features - 1);
         $this_count >= 0;
         --$this_count) {
    my $this_entry_feature = $this_entry_overlap_features[$this_count];
  
    for (my $next_count = (@next_entry_overlap_features - 1);
         $next_count >= 0;
         --$next_count) {
      my $next_entry_feature = $next_entry_overlap_features[$next_count];
  
      # check each /gene in $this_entry_feature against each /gene in
      # $next_entry_feature.  any match means that the features should
      # be merged
  
      my @this_feature_gene_names =
        map {
          $_->value
        } grep {
          $_->name eq "gene"
        } $this_entry_feature->qualifiers;
  
      my @next_feature_gene_names =
        map {
          $_->value
        } grep {
          $_->name eq "gene"
        } $next_entry_feature->qualifiers;
  
      for my $this_feature_gene_name (@this_feature_gene_names) {
  
        if (grep {lc $_ eq lc $this_feature_gene_name}
            @next_feature_gene_names) {
  
          ##########
  
          #  merge the two features by removing $next_entry_feature and
          #  adding stuff to $this_entry_feature
  
          # merge the locations
  
          my $new_location =
            merge_locations ($this_entry_feature->location,
                             $next_entry_feature->location);
  
          if (! defined $new_location) {
            # locations can't be merged probably means that one of the
            # sequences is in the wrong orientation
            last this_entry_loop;
          }
  
          # remove $this_entry_feature from @this_entry_overlap_features
  
          splice @this_entry_overlap_features, $this_count, 1;
  
          # remove $next_entry_feature from @next_entry_overlap_features
  
          splice @next_entry_overlap_features, $next_count, 1;
  
          $this_entry_feature->location ($new_location);
  
          # merge the qualifiers from the $next_entry_feature into
          # $this_entry_feature
  
          add_qualifiers ($next_entry_feature, $this_entry_feature);
  
          # remove /partial and /codon_start qualifiers
          for (my $i = scalar($this_entry_feature->qualifiers) - 1;
               $i >= 0;
               --$i) {
            my @qualifiers = @{$this_entry_feature->qualifierList};
            if ($qualifiers[$i]->name =~ /^(partial|codon_start)\s*$/) {
              splice @{$this_entry_feature->qualifierList}, $i, 1;
            }
          }
  
          # add the feature to $entry_1
          $ft_1->addFeature ($this_entry_feature);
  
          next this_entry_loop;
        }
      }
    }
  }
  
  
  # add the remaining features from entry 2 to entry 1
  
  for my $feature ($ft_2->features) {
    $ft_1->addFeature ($feature, "end");
  }
  
  
  # add the left over overlap CDSs
  
  for my $feature (@this_entry_overlap_features, 
                   @next_entry_overlap_features) {
    $ft_1->addFeature ($feature, "end");
  }
  
  
  #############
  
  #  now show /partial CDS features
  
  my $partial_count = 0;
  
  for my $feature ($ft_1->features) {
    if (grep {$_->name =~ /^partial\s*$/} @{$feature->qualifierList}) {
      ++$partial_count;
    }
  }
  
  print STDERR <<EOF;
after joining there are $partial_count partial CDS features left
EOF
}


my $first_file = shift;

print STDERR "first entry: $first_file\n";

open IN_FILE_1, "$first_file" or die "$0: can't open $first_file\n";

my $set_1 = Read Embl(\*IN_FILE_1);

my @entries_1 = $set_1->entries();

if (@entries_1 != 1) {
  die "$0: $ARGV[0] does not contain one entry\n";
}

my $entry_1 = pop @entries_1;


for my $this_file (@ARGV) {
  print STDERR "attempting to add $this_file to the contig\n";

  open IN_FILE_2, "$this_file" or die "$0: can't open $this_file\n";
  
  my $set_2 = Read Embl(\*IN_FILE_2);

  my @entries_2 = $set_2->entries();

  if (@entries_2 != 1) {
    die "$0: $ARGV[1] does not contain one entry\n";
  }

  my $entry_2 = pop @entries_2;

  open OVERLAP_TAB, ">overlap.tab" or
    die "cannot open overlap.tab for writing\n";
  
  merge_entries ($entry_1, $entry_2);

  close OVERLAP_TAB;
}

# now write the entry
  
$entry_1->toStream(\*STDOUT);
