# -*-Perl-*- 
# 
# $Id: Assembly.pm,v 1.17 1998/09/28 10:04:55 dn1 Exp $
#
# $Log: Assembly.pm,v $
# Revision 1.17  1998/09/28  10:04:55  dn1
# warn woah becomes optional
#
# Revision 1.16  1998/02/04  17:13:29  sd
# Removed next after 'Woah' so that something useful is gained
# from problem reads.
#
# Revision 1.15  1998/02/04  12:11:59  sd
# Oops. examineRead now returns non-zero.
#
# Revision 1.14  1998/02/02  14:12:57  sd
# No longer use ReadTypes->isTerminator
#
# Revision 1.13  1998/01/08  14:56:23  sd
# Sensible changes
#
# Revision 1.12  1997/07/29  15:53:44  dn1
# Stolen reads can be coped with now. Unresolved consenus regions
# fixed (was out by 1).
#
# Revision 1.11  1997/07/21  09:44:59  apm
# more Perl5.004 fixing!
#
# Revision 1.10  1997/07/21  09:01:37  apm
# Fixed Perl5.004 issues
#
# Revision 1.9  1997/06/24  11:31:00  dn1
# new version
#
# Revision 1.8  1997/04/24  08:53:44  dn1
# Added copyright
#
# Revision 1.7  1996/07/08  08:27:03  dmh
# Supports extraction of consensus base quality values.
#
# Revision 1.6  1996/06/25  09:18:13  dmh
# Added getConsensusTags
#
# Revision 1.5  1996/04/04  16:11:12  dmh
# Added tag coverage croutines
#
# Revision 1.1.1.1  1996/03/27  14:40:01  dmh
# Fixed for new usage of ReadTypes.pm
#
# Revision 1.4  1996/03/26  12:12:23  dmh
# More minor fiddles.
#
# Revision 1.2  1996/03/19  18:20:57  dmh
# Uses Cover.pm
#
# Revision 1.1  1996/03/18  13:55:44  badger
# Initial revision
#
#
# 
#######################################################################
# This software has been created by Genome Research Limited (GRL).    # 
# GRL hereby grants permission to use, copy, modify and distribute    # 
# this software and its documentation for non-commercial purposes     # 
# without fee at the user's own risk on the basis set out below.      #
# GRL neither undertakes nor accepts any duty whether contractual or  # 
# otherwise in connection with the software, its use or the use of    # 
# any derivative, and makes no representations or warranties, express #
# or implied, concerning the software, its suitability, fitness for   #
# a particular purpose or non-infringement.                           #
# In no event shall the authors of the software or GRL be responsible # 
# or liable for any loss or damage whatsoever arising in any way      # 
# directly or indirectly out of the use of this software or its       # 
# derivatives, even if advised of the possibility of such damage.     #
# Our software can be freely distributed under the conditions set out # 
# above, and must contain this copyright notice.                      #
#######################################################################


package Assembly;

sub setnowoah {
    ($nowoah) = @_;
    #warn "nowoah -> $nowoah \n";
}


require AceTools;
require Cover;
require ReadTypes;

#
# Global object data struct
#
$contigs = {};


##########################################
#
# Construct....
#
##########################################

sub new {

    my $type = shift;
    my $self = {@_};

    bless $self;

}

##########################################
#
# Initialise....(read in the ace file)
#
##########################################
sub init {

    my $type = shift;

    $type->{'ace'} = AceTools->new($type->{'cafile'});

    my $contig;

    foreach $contig ($type->{'ace'}->theContigs) {
	$contigs->{$contig} = {} if (! defined $contigs->{$contig});
	$contigs->{$contig}->{'type'}='contig';
    }

}


#########################################################
#
# Main method section
#
#########################################################

sub getContigList {

    return (keys(%$contigs));

}


sub getContigLength {

    my ($type, $contig)=@_;

    # 
    # If length for this contig not defined, read
    # them all. This should only happen once.
    # 
    unless (defined($contigs->{$contig}->{'length'})) {
	my $hash_ref = $type->{'ace'}->readContigLengths();	     
	foreach (keys(%$contigs)) {
	    $contigs->{$_}->{'length'} = $$hash_ref{$_};
	}
    }
    
    # 
    # Lengths should now all be defined in $contig
    # 
    unless (defined($contigs->{$contig}->{'length'})) {
	warn "Couldn't get length for $contig\n";
    } else {
	#print "Contig: $contig Length: ", $contigs->{$contig}->{'length'}, "\n";
	return $contigs->{$contig}->{'length'};
    }

}

sub getAfEntries {

    my ($type, $contig)=@_;

    unless (defined $contigs->{$contig}->{'af'}) {
	$contigs->{$contig}->{'af'} = $type->{'ace'}->getAfEntries($contig);
    }
    return $contigs->{$contig}->{'af'};

}

sub getReadsPerContig {

    my ($type, $contig)=@_;

    my $af_ref = $type->getAfEntries($contig);

    my @reads;
    my $i;
    for($i=0 ; $i<=$#{$af_ref} ; $i++) {
	push (@reads, $$af_ref[$i]->[0]);
    }

    return (@reads);

}


sub getTotalContigLength {

    my ($type)=@_;
        
    my $total = 0;
    foreach (keys(%$contigs)) {
	$total += $type->getContigLength($_);
    }

    $total;

}


sub get_strand_info {
    my ($type, $contig, $info)=@_;

     unless (defined($contigs->{$contig}->{$info})) {
	 $type->defineStrandInfo($contig);
    }

    return $contigs->{$contig}->{$info};
}



sub defineStrandInfo {

    my ($type, $contig)=@_;

    my $af_ref = $type->getAfEntries($contig);

    my $primer_fwd;
    my $primer_rev;
    my $term_fwd;
    my $term_rev;

    $contig_len = $type->getContigLength($contig);

    $term_fwd = Cover->new(1, $contig_len);
    $term_rev = Cover->new(1, $contig_len);
    $primer_fwd = Cover->new(1, $contig_len);
    $primer_rev = Cover->new(1, $contig_len);

    my $i;
    for($i=0 ; $i<=$#{$af_ref} ; $i++) {
	    
	my ($read, $c1,$c2) = ($$af_ref[$i]->[0], $$af_ref[$i]->[1], $$af_ref[$i]->[2]);

	########################################
	my $read_info=ReadTypes->new($read);
	if (! $read_info->get('good')) {
	    if ($nowoah == 0){
		warn "Woah! check read $read\n";
	    }
	}
	
	########################################

	if ($read_info->getItemInfo('dye') eq "terminator") {
	    if ( $c1 < $c2 ) {
	      $term_fwd->cover($c1, $c2);
	    } else {
	      $term_rev->cover($c1, $c2);
	    }
	}  else {
	    if ( $c1 < $c2 ) {
	      $primer_fwd->cover($c1, $c2);
	    } else {
	      $primer_rev->cover($c1, $c2);
	    }
	}
	
    }

    $term_fwd->complement();
    $term_rev->complement();
    $primer_fwd->complement();
    $primer_rev->complement();

    #########################################################
    
    my $term_fwd_only = $term_fwd->copy()->exclude($term_rev)->exclude($primer_fwd)->exclude($primer_rev);
    my $term_rev_only = $term_rev->copy()->exclude($term_fwd)->exclude($primer_fwd)->exclude($primer_rev);
    my $primer_fwd_only = $primer_fwd->copy()->exclude($term_fwd)->exclude($term_rev)->exclude($primer_rev);
    my $primer_rev_only = $primer_rev->copy()->exclude($term_fwd)->exclude($term_rev)->exclude($primer_fwd);
    
    


    $contigs->{$contig}->{'primer_fwd'}      = $primer_fwd;
    $contigs->{$contig}->{'primer_rev'}      = $primer_rev;
    $contigs->{$contig}->{'term_fwd'}        = $term_fwd;
    $contigs->{$contig}->{'term_rev'}        = $term_rev;
    $contigs->{$contig}->{'primer_fwd_only'} = $primer_fwd_only;
    $contigs->{$contig}->{'primer_rev_only'} = $primer_rev_only;
    $contigs->{$contig}->{'term_fwd_only'}   = $term_fwd_only;
    $contigs->{$contig}->{'term_rev_only'}   = $term_rev_only;

}


#########################################################
#
# Draw clumps of a tag type.
# 
#########################################################
sub getTagRegions {

    my ($type, $contig, $tagtype)=@_;

    ##################################
    # 
    # Get contig length and initialise cover array
    # 
    ##################################

    my $contig_len = $type->getContigLength($contig);

    my $tags;
    $tags = Cover->new(1, $contig_len);

    ##################################
    # 
    # Read tags from contig object
    # 
    ##################################

    my $af_ref = $type->getAfEntries($contig);

    my $i;
    for($i=0 ; $i<=$#{$af_ref} ; $i++) { # Read each AF line
	
	my ($read, $c1,$c2, $r1, $r2) = ($$af_ref[$i]->[0], 
					 $$af_ref[$i]->[1], 
					 $$af_ref[$i]->[2],
					 $$af_ref[$i]->[3], 
					 $$af_ref[$i]->[4]);

	
	##################################
	# 
	# Get a table of tag coords for requested type
	# 
	##################################

	my $taglist_ref;
	
	$taglist_ref = $type->{'ace'}->getTagCoords($read, $tagtype);
	if ($tagtype eq "STOL") { #
	    if ($type->{'ace'}->getSequence($read) =~ /\nStolen/) {
		push(@$taglist_ref,[ $c1, $c2 ]);
	    }
	}	
	##################################
	# 
	# Cover tagged areas
	# 
	##################################

	foreach (@$taglist_ref) {
	    my ($t1,$t2)=@$_;
	    if ($c1<$c2) {
		($tc1, $tc2) = ($c1+($t1-$r1), $c2+($t2-$r2));
	    } else {
		($tc1, $tc2) = ($c1-($t1-$r1), $c2-($t2-$r2));
	    }

	    ($tc1, $tc2) = ($tc1 < $tc2) ? ($tc1, $tc2) : ($tc2, $tc1);
	    
 	    $tc1 = ($tc1<1) ? 1 : $tc1;
 	    $tc2 = ($tc2<1) ? 1 : $tc2;
 
 	    $tc1 = ($tc1>$contig_len) ? $contig_len : $tc1;
 	    $tc2 = ($tc2>$contig_len) ? $contig_len : $tc2;

# 	    print "$read $c1 $c2 $r1 $r2   $t1 $t2\n";
# 	    print "COVER: from: $tc1, to: $tc2\n";

	    $tags->cover($tc1, $tc2);
	}
    }

    $tags->complement();

    return $tags;

}


#########################################################
#
# Get list of consensus tags straight
# from caf file.
# 
#########################################################
sub getConsensusTags {

    my ($this, $contig, $tagtype)=@_;

    return $this->{'ace'}->getConsensusTags($contig, $tagtype);

}


#########################################################
#
# Get confidence values for a contig.
# 
#########################################################
sub getConsensusQuality {

    my ($this, $contig)=@_;

    return $this->{'ace'}->getConsensusQuality($contig);

}

#########################################################
#
# Return ref to list of regions in consensus
# containing N's.
# 
#########################################################

sub getUnresolved {

    my ($type, $contig)=@_;

    my $seq = $type->{'ace'}->getDNA($contig);

    my $i;
    my $start=-1;
    my @nlist;

    for ($i=0 ; $i < length($seq) ; $i++) {
	my $base = substr($seq, $i, 1);

	if ($base =~ /[nN\-]/) {
	    $start = $i if $start == -1;
	} elsif (-1 != $start) {
	    #print "Contig: $contig N: $start, $i\n";
	    push(@nlist, [$start+1, $i]);
	    $start = -1;
	}
    }
    if (-1 != $start) {
	#print "Contig: $contig N: $start, $i\n";
	push(@nlist, [$start+1, $i]);
    }

    return \@nlist;

}



############################################################
#
# Get template coverage 
# 
############################################################

sub getTemplateCoverage {

    my ($type, $contig)=@_;

    my $af_ref = $type->getAfEntries($contig);
    my @ordered_af = $type->bothEndsAf($af_ref);
    my @singles = $type->findSingles(\@ordered_af);

    return (@singles);

}


########################################################
# bothEndsAf
# 
# Sort Af table into positional order. Each line appears
# twice - start points and end points.
#
########################################################
sub bothEndsAf {

    my ($type, $af_ref) = @_;

    my %mark = ();
    my @order = ();
    my $high = 0;


    foreach $af (@$af_ref) {

	my ($read, $c1, $c2) = @$af;
	
	if ($c1 > $high) {
	    push(@{$mark{$c1}}, [$read, $c1, $c2]);
	}

	if ($c2 > $high) {
	    push(@{$mark{$c2}}, [$read, $c1, $c2]);
	}
    }

    foreach $index (sort{$a<=>$b}(keys(%mark))) {
	foreach (@{$mark{$index}}) {
	    ($read, $c1, $c2) = @$_;
	    push(@order, [ $read, $c1, $c2 ]);
	}
    }
    
    undef(%mark);

    return (@order);
    
}


sub findSingles {

    my ($type, $order_ref) = @_;

    my $pos = 0;
    my $state = 0;
    my @stack = ();
    my ($start, $end);
    my @single = ();

    #warn "\n\n" if $debug;

    foreach (@$order_ref) {

	my ($read, $c1, $c2) = @$_;
	#warn "NEXT: $read, $c1, $c2\n" if $debug;

	my ($lo, $hi) = sort {$a<=>$b}  ($c1, $c2);

	unless (grep($_->[0] eq $read, @stack)) {
 	    #warn "PUSH: $read, $lo $hi\n" if $debug;
	    push(@stack, [$read, $lo, $hi]);
	    $pos = $lo;
	} else {
	    $pos = $hi;
	}
	
	#warn "POS: $pos\n" if $debug;

	my %seen;
	my $poped = 0;
	for ($j=0 ; $j <= $#stack ; $j++) {	# Look thru stack
	    
	    #warn "STACK COUNT: ", $#stack, "COUNT: $j\n" if $debug;

	    my ($s_read, $from, $to) = ($stack[$j]->[0], $stack[$j]->[1], $stack[$j]->[2]);
	    
	    #warn "STK: $s_read, $from, $to\n" if $debug;
	    
	    if ($to == $pos && $poped == 0) {
		#warn "At pos: $pos POP: $s_read, $from, $to\n" if $debug;
		splice(@stack, $j, 1);
		$j--;
		$poped =1;
		next;
	    } 

	    my $read_info=ReadTypes->new;
	    $read_info->examineRead($s_read);
	    my $template=$read_info->getItemInfo('template');
	    if ($template eq '') {
		#warn "WARNING: Couldn't derive template name for read: \'$s_read\'\n";
	    }
	    #warn "TEMPLATE: $template\n" if $debug;
	    $seen{$template}++;
	    
	}
	
	my $t_count = scalar(keys(%seen));
 	#warn "TEMPS: $t_count STATE: $state\n" if $debug;
	
	if ($t_count <= 1 && $state == 0) { # Only one template
	    $state=1;
	    if ($pos==$lo){ $start=$pos; }
	    else { $start=$pos+1; }
 	    #warn "START: $start\n" if $debug;
	} elsif ($t_count > 1 && $state == 1) {
	    $state=0;
	    $end=$pos-1;
 	    #warn  "END: $end\n" if $debug;
	    push(@single, [$start, $end]) if ($start<=$end);
	}
	
    }

    if ( $state == 1) {
	$end=$pos;
 	#warn "END: $end\n" if $debug;
	push(@single, [$start, $end]) if ($start<=$end);
    }

    return (@single);

}



sub min {

    my $minval = shift;

    my $cmp;
    while($cmp = shift) {
	if ($minval > $cmp) { $minval = $cmp; }
    }

    $minval;

}


sub byReadOrder {
    min($a->[1],$a->[2]) <=> min($b->[1],$b->[2]);
}

sub printAf {

    my ($type, $af_ref)=@_;

    my $i;
    for($i=0 ; $i<=$#{$af_ref} ; $i++) {
	print $$af_ref[$i]->[0], $$af_ref[$i]->[1], $$af_ref[$i]->[2], "\n";
    }

}


1;














