#!/usr/bin/perl -w
# -*-CPerl-*-
# 	$Id: pmmulti.pl,v 1.6 2003/12/22 08:42:12 ivo Exp $	

use strict;
use warnings;
use Pod::Usage;
use Getopt::Long qw(:config pass_through);
use vars qw($VERSION);
use FindBin qw($Bin);

$VERSION = 1.1;

my $fast = 0;
my $slow = 0;
my $man = 0;

sub usage {
  die "Usage: $0 [--fast_progressive] [--slow_pairwise] [pmcomp options]\n";
}

my $pairopt = '';
my @ToDo;

sub pairopt { # pass on argument to pmcomp.pl
  $pairopt .= ' -' . join(' ', @_);
}
GetOptions("fast_progressive!" => \$fast,
	   "slow_pairwise!" => \$slow,
	   "man" => \$man,
	   "g=f" => \&pairopt,
	   "D=i" => \&pairopt,
	   "seqw=f" => \&pairopt,
	   "nolog" => \&pairopt,
	   "fast!" => \&pairopt,
	   "endgaps!" => \&pairopt,
	   '<>' => sub {push @ToDo, @_}
	  ) or pod2usage(2);

pod2usage(-exitstatus => 0, -verbose => 2) if $man;


my $PMCOMP = "$Bin/pmcomp.pl";  # search pmcomp in same directory as pmmulti
my @PairAlignments;
my $WORKDIR = './pm_tmp/';
my @SM;
`rm -rf $WORKDIR`;  # Just in case and to help debugging, remove everything
`mkdir $WORKDIR`;

print "# Input files\n";

if (!@ToDo) {
  while (<>) { # alternatively read input names from stdin
    chomp;
    push  @ToDo, $_;
  }
}
my $NN=0;
foreach (@ToDo) {
  s/_dp.ps$//;
  $NN++;
  my $f = $_;
  $f .= '_dp.ps' unless /.ps$/;
  printf("# %3d <- %s\n",$NN, $f);
  `cp $f $WORKDIR/$NN.ps`;
}

print "# BE PATIENT ... computing pairwise alignments ";
print "using ", $slow ? "slow/accurate" : "fast/profile", " method\n";

# First compute all pairwise alignments, the scores go into @SM;

$pairopt .= $slow ? " --nofast" : " --fast";
for (my $i=1; $i<=$NN; $i++) {
  for (my $j=$i+1; $j<=$NN; $j++) {
    print "#  $i with $j : score = ";
    my $fna =  "$WORKDIR" . "/X" . "$i" ."_" . "$j" . ".aln" ;
    my $fnp =  "$WORKDIR" . "/X" . "$i" ."_" . "$j" . "_aln_dp.ps" ;
    `$PMCOMP @ARGV $pairopt $WORKDIR/$i.ps   $WORKDIR/$j.ps  > $fna`;
    `mv aln_dp.ps $fnp `;
    open(AL, "<$fna")  || die "can't read $fna $!";
    <AL>; $_ = <AL>;
    warn "can't parse score" unless /(-?\d+\.?\d*)/;
    $SM[$i][$j] = $1;
    close(AL);
    $SM[$j][$i] = $SM[$i][$j];
    print $SM[$i][$j] ."\n";
  }
}

print "#\n# WPGMA guide tree and progressive multiple alignment\n";

# Compute profile alignments

my $profopt = $fast ? "--fast" : "--nofast";
my (@und, $im, $jm, $k, $xmin);
for (my $i=1; $i<=$NN; $i++) { $und[$i]=1; }
for (my $i=$NN+1; $i<=2*$NN; $i++) { $und[$i]= (-1);}

for ($k=$NN+1; $k<2*$NN; $k++) {
  $xmin = 0.0;
  for (my $i=1; $i<=2*$NN; $i++) {
    next if $und[$i]<0;
    for (my $j=1; $j<=2*$NN; $j++) {
      next if $und[$j]<0;
      next if $i == $j;
      if ($xmin < $SM[$i][$j] ) { $xmin =  $SM[$i][$j]; $im = $i; $jm = $j; }
    }
  }
  print "# Merge $im ($und[$im]) with $jm ($und[$jm]) -> $k at $xmin ...";
  my $fna =
    "$WORKDIR" . "/X" . "$im" ."_" . "$jm" . "=" . "$k" . ".aln" ;
  my $fnp =
    "$WORKDIR" . "/X" . "$im" ."_" . "$jm" . "=" . "$k" . "_aln_dp.ps" ;
  `$PMCOMP $profopt $WORKDIR/$im.ps   $WORKDIR/$jm.ps  > $fna`;
  `mv aln_dp.ps  $WORKDIR/$k.ps `;
  open(AL,"<$fna") || die "can't read from $fna $!";
  <AL>; $_ = <AL>;
  warn "can't parse score" unless /(-?\d+\.?\d*)/;
  my $sco = $1;
  my @ALI;
  while (<AL>) {
    chomp;
    push @ALI, $_ if /^[AUGCRYN(.)-]+/;
  }
  close(AL);
  print "score = $sco\n";

#  foreach (@ALI) { print "$_\n"; }

  push @PairAlignments, [$k, $im, $jm, @ALI];

  for (my $l=1; $l<=2*$NN; $l++) {
    next if $und[$l]<0;
    next if $l == $k;
    next if $l == $im;
    next if $l == $jm;
    if (defined($SM[$l][$im]) && defined($SM[$l][$jm])) {
      $SM[$k][$l] =
	($SM[$l][$im])*($und[$im]/($und[$im]+$und[$jm])) +
	  ($SM[$l][$jm])*($und[$jm]/($und[$im]+$und[$jm]));
    }
    else { print "undef: $l, $im, $jm \n" ;}
    $SM[$l][$k] = $SM[$k][$l];
  }
  $und[$k] =$und[$im]+$und[$jm];
  $und[$im]= (-1);
  $und[$jm]= (-1);
}

#foreach (@PairAlignments) { print "$_\n"; }

my $merges = @PairAlignments;

# print $merges . "\n";
my $m = $merges-1;
# print "\n:$m\n";
my @MultipleAlignment;

print "# Finally we re-assemble the multiple alignment from the pairwise
ones \n#\n";

push @MultipleAlignment, @{$PairAlignments[-1]}[1,3,4,2,5,6];


# foreach (@MultipleAlignment) { print "$_\n"; } print "\n";

for ($m=$merges-2; $m>=0; $m--) {
  my ($which, $which1, $which2) = @{$PairAlignments[$m]}[0..2];

  for (my $i=0; $i<scalar(@MultipleAlignment); $i++) {
    if ($MultipleAlignment[$i] eq $which) {
      my $ali  = $MultipleAlignment[$i+1];
      my ($id1, $id2, $seq1, $str1, $seq2, $str2) =
	@{$PairAlignments[$m]}[1..6];

      # now use $ali produce the additional gaps in $str*, $seq*
      my ($seq1a, $seq2a, $str1a, $str2a) = ('','','','');

      for (my $j=0, my $j1=0 ; $j<length($ali); $j++) {
	if ( substr($ali,$j,1) eq '-') {
	  $seq1a = $seq1a . "-" ;
	  $seq2a = $seq2a . "-" ;
	  $str1a = $str1a . "-" ;
	  $str2a = $str2a . "-" ;
	}
	else {
	  $seq1a = $seq1a . substr($seq1, $j1, 1);
	  $seq2a = $seq2a . substr($seq2, $j1, 1);
	  $str1a = $str1a . substr($str1, $j1, 1);
	  $str2a = $str2a . substr($str2, $j1, 1);
	  $j1++;
	}
      }

      splice @MultipleAlignment, $i, 3,
	$id1, $seq1a, $str1a, $id2, $seq2a, $str2a;

      last;
    }
  }
#  print "$m:\n";
#  foreach (@MultipleAlignment) { print "$_\n"; }

}

for (my $i=0; $i<=$merges; $i++) {
  printf("%-16s  %s\n", $ToDo[$MultipleAlignment[$i*3]-1],
	 $MultipleAlignment[$i*3+1]);
  printf("    %-12d  %s\n", $MultipleAlignment[$i*3],
	 $MultipleAlignment[$i*3+2]);
}

my $cseq = '';
my $cstr = '';
my ($seq1, $str1, $seq2, $str2) = @{$PairAlignments[-1]}[3..6];
# most of this unused

for (my $i=0; $i<length($seq1); $i++) {
  my $c1 = substr($str1,$i,1);
  my $c2 = substr($str2,$i,1);

  if ( $c1 eq $c2 ) { $cstr = $cstr . $c1; next; }
  if ( $c1 eq '-')  { $cstr = $cstr . '~'; next; }
  if ( $c2 eq '-')  { $cstr = $cstr . '~'; next; }
  $cstr = $cstr . '?'; # do something as default
}

for (my $i=0; $i<length($seq1); $i++) {
#
# Replace by a direct computation of consensus from sequences
#  lateron
#
  my $n1 = substr($seq1,$i,1);
  my $n2 = substr($seq2,$i,1);
  if ( $n1 eq $n2 ) { $cseq = $cseq . $n1; next; }
  if ( $n1 eq '-')  { $cseq = $cseq . $n2; next; }
  if ( $n2 eq '-')  { $cseq = $cseq . $n1; next; }
  $cseq = $cseq . 'N'; # do something as default
}

printf("CONSENSUS_SEQ     %s\n", $cseq);
printf("CONSENSUS_STR     %s\n", $cstr);

__END__

=head1 NAME

pmmulti.pl - Multiple alignment of RNA pair matrices

=head1 SYNOPSIS

pmmulti.pl [pmmulti opts] [pmcomp opts]

Options:
        --fast         fast/approximate progressive alignments
        --slow         slow/accurate pairwiser alignments
        --help         brief help message
        --man          full documentation

=head1 OPTIONS

=over 8

=item B<--fast_progressive>

use the fast but inaccurate profile alignment for the progressive
multiple alignments.  Default is to use the slow (Sankhoff) method.

=item B<--slow_pairwise>

use Sankhoff mehtod even for the pairwise alignments.  Default is
to use the fast string-like profile alignment for this stage.

=item B<other options>

any unrecognized options will be passed on to pmcomp.pl

=back

=head1 DESCRIPTION

B<pmmulti.pl> computes mutliple alignments of RNA pair probability
matrices using a variant of the Sankhoff algorithm. It is implemented
as a wrapper that calls C<pmcomp.pl> to actually compute pairwise
alignments.

Input consists simply of a list of sequence names, the current
directory should contain a dot plot files with pair probabilities for
each sequence named in the input in the foprmat produced by C<RNAfold
-p>.

A directory C<pm_tmp> in the current directory will be created to
store intermediate results. Any file in an existing C<pm_tmp>
directory may be overwritten.

=cut
