#!/usr/local/bin/perl
# -*-Perl-*- 
#######################################################################
# 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.                      #
#######################################################################
# 
# $Id: Caftools.pm,v 1.13 2001/06/29 11:12:56 mng Exp $

package Caftools;

use Carp;
use strict;
use ReadTypes;
use Uncovered;
use Exptools;
use FileHandle;
use TempFiles;

my $tFiles = TempFiles->new(prefix => 'caftools');

sub new {
    my $type = shift;
    my %params = @_;
    my $self = {};
    $self->{verbose} = $params{verbose};
    $self->{CAF} = {DNA => {}, BaseQuality => {}};
    return bless $self, $type;
}

sub readCafFile {
    my ($self, $file, $duplicatePolicy) = @_;
    confess "Wrong type" unless ref $self;
    $duplicatePolicy = 'old' unless $duplicatePolicy;
    my $fileRef = (ref($file)) ? $file : (FileHandle->new("< $file") || die "Couldn't open $file for reading\n");
    my $size = -s $fileRef;
    $size = 1 unless ($size);
    my ($name, $class);
    my ($oldPercent, $verbose) = (-1, $self->{verbose});
    print STDERR "Reading CAF file\n" if ($verbose);
    my @x;
    while (<$fileRef>) {
	if ($verbose) {
	    my $percent = int(tell($fileRef) * 100 / $size);
	    print STDERR "\rRead $percent\%" if ($percent > $oldPercent);
	    $oldPercent = $percent;
	}
        s/\/\/.*// ;			# ace file comments
        if ( /^DNA\s+(:\s+)?(\S+)/ ) {
            $name = $2;
	    @x = ();
	    while (<$fileRef>) {
		last if /^\s*$/ ;
		s/\s+//g ;		# remove whitespace
		push(@x,$_);
	    }
	    $self->{CAF}->{DNA}->{$name} = join("",@x)
		unless ($self->{CAF}->{DNA}->{$name} && $self->_reportDuplicate($name, 'DNA', $duplicatePolicy));
	} elsif ( /BaseQuality\s+(:\s+)?(\S+)/ ) {
            $name = $2;
	    @x = ();
	    while (<$fileRef>) {
		last if /^\s*$/;
		chomp;
		push(@x, $_);
	    }
	    $self->{CAF}->{BaseQuality}->{$name} = join(" ", @x)
		unless ($self->{CAF}->{BaseQuality}->{$name} && $self->_reportDuplicate($name, 'base quality', $duplicatePolicy));
	} elsif ( /^(\S+)\s+(:\s+)?(\S+)/ ) {
	    $class = $1;
	    $name = $3;
	    $self->{CAF}->{$class} = {} unless $self->{CAF}->{$class};
	    @x = ();
	    while (<$fileRef>) {
		last if /^\s*$/ ;
		s/\\\n$// ;		# continuation lines
		push(@x,$_);
	    }
	    if (defined($self->{CAF}->{$class}->{$name})) {
	        $self->{CAF}->{$class}->{$name} .= join("",@x) unless $self->_reportDuplicate($name, $class, $duplicatePolicy);
	    } else {
	        $self->{CAF}->{$class}->{$name} = join("",@x);
	    }
	}
	last if eof($fileRef);
    }
    print STDERR "\rRead 100%\n" if ($verbose);
}

sub _reportDuplicate {
    #return 1 if duplicate is to be ignored, 0 if not (does whatever old readAceFile did).
    my ($self, $name, $class, $duplicatePolicy) = @_;
    confess "Wrong type" unless ref($self);
    if ($duplicatePolicy eq 'die') {
	die "Fatal error: Found duplicate $class entry $name in CAF file.\n";
    } elsif ($duplicatePolicy eq 'warn') {
	warn "Found duplicate $class entry $name in CAF file.\n";
	return 1;
    } elsif ($duplicatePolicy eq 'ignore') {
	return 1;
    } elsif ($duplicatePolicy eq 'old') {
	return 0;
    } else {
	croak "Bad policy for handling CAF file duplicate reads $duplicatePolicy";
    }
}

sub writeCafFile {
    my ($self, $file) = @_;
    confess "Wrong type" unless ref($self);
    my $verbose = $self->{verbose};
    my $fileRef = (ref($file)) ? $file : (FileHandle->new("> $file") || die "Couldn't open $file for writing\n");
    my $className;
    foreach $className (keys %{$self->{CAF}}) {
	my $class = $self->{CAF}->{$className};
	print STDERR "Writing class $className\n" if ($verbose);
	my @classKeys = keys %{$class};
	my ($size, $written, $oldPercent) = (scalar(@classKeys), 0, -1);
	my $name;
	foreach $name (@classKeys) {
	    if ($verbose) {
		my $percent = int($written++ * 100 / $size);
		print STDERR "\rWritten $percent\%" if ($percent > $oldPercent);
		$oldPercent = $percent;
	    }
	    next unless ($class->{$name});
	    print $fileRef "\n$className : $name\n" || die "Error in writeCafFile while writing CAF file\n";
	    if ($className eq "DNA") {
		my $DNA = $class->{$name};
		my $len = length($DNA);
		$DNA =~ s/[^ACGTUacgtu\-\*]/N/g; # translate all non-nucleotides to N's
		my $i;
		for ($i = 0; $i < $len; $i += 60) {
		    my $line = substr($DNA, $i, 60);
		    print $fileRef "$line\n" || die "Error in writeCafFile while writing CAF file\n";
		} 
	    } else {
		print $fileRef "$class->{$name}\n" || die "Error in writeCafFile while writing CAF file\n";
	    }
	}
	print STDERR "\rWritten 100\%\n" if ($verbose);
    }
    print $fileRef "\n" || die "Error in writeCafFile while writing CAF file\n";
}

sub getEntry {
    my ($self, $entryName) = @_;
    my $entry = {};
    my $class;
    foreach $class (keys %{$self->{CAF}}) {
	$entry->{$class} = $self->{CAF}->{$class}->{$entryName};
    }
    return $entry;
}

sub get_entry_ref {
    my ($self, $entry_name) = @_;

    my %entry;
    while (my ($class, $entries) = each(%{$self->{CAF}})) {
	$entry{$class} = \$entries->{$entry_name};
    }
    return \%entry;
}

sub insertEntry {
    my ($self, $entry, $entryName) = @_;
    my $class;
    foreach $class (keys %{$entry}) {
	$self->{CAF}->{$class} = {} unless ($self->{CAF}->{$class});
	$self->{CAF}->{$class}->{$entryName} = $entry->{$class};
    }
}

sub entriesList {
    my ($self) = @_;
    confess "Wrong type" unless ref($self);
    return keys(%{$self->{CAF}->{Sequence}});
}

sub contigsList {
    my ($self) = @_;
    confess "Wrong type" unless ref($self);
    return grep($self->{CAF}->{Sequence}->{$_} =~ /Is_contig\n/, keys(%{$self->{CAF}->{Sequence}}));
}

sub readsList {
    my ($self) = @_;
    confess "Wrong type" unless ref($self);
    return grep($self->{CAF}->{Sequence}->{$_} =~ /Is_read\n/, keys(%{$self->{CAF}->{Sequence}}));
}

sub assembliesList {
    my ($self) = @_;
    confess "Wrong type" unless ref($self);
    return grep($self->{CAF}->{Sequence}->{$_} =~ /Is_assembly\n/, keys(%{$self->{CAF}->{Sequence}}));    
}

sub trashRead {
    my ($self, $read) = @_;
    confess "Wrong type" unless ref($self);
    delete $self->{CAF}->{Sequence}->{$read};
    delete $self->{CAF}->{DNA}->{$read};
    delete $self->{CAF}->{BaseQuality}->{$read};
}


sub trashContig {
  my($self,$contig) = @_;
  confess "Wrong type" unless ref($self);
  
  foreach my $line (split("\n", $self->{CAF}->{Sequence}->{$contig})) {    
    my ($af, $read) = (split(/\s+/, $line))[0,1];
    next unless ($af eq "Assembled_from");    
    $self->trashRead($read);
  }
  $self->trashRead($contig);
}

sub trashAll {
    my ($self) = @_;
    confess "Wrong type" unless ref($self);
    $self->{CAF} = {DNA => {}, BaseQuality => {}};
}

sub tag {
    my ($self, $contig, $left, $right, $newTag, $comment, $exclude, $readExtentsCache) = @_;
    confess "Wrong type" unless ref($self);
    warn "Tagging reads for $newTag in contig $contig from $left through $right\n" if $self->{debug};
    my $leftEnd = ($left <= $right)? $left : $right;
    my $rightEnd = ($left <= $right)? $right : $left;
    my $tagKey = ($newTag eq "CVEC") ? "Clone_vec":"Tag";

    # Tag consensus
    $self->{CAF}->{Sequence}->{$contig} .= "$tagKey $newTag $left $right \"$comment\"\n";

    # Tag readings (more difficult)
    my $readExtents;
    if ($readExtentsCache) {
	unless ($readExtentsCache->{$contig}) {$readExtentsCache->{$contig} = $self->_getReadExtents($contig);}
	$readExtents = $readExtentsCache->{$contig};
    } else {
	$readExtents = $self->_getReadExtents($contig);
    }
    my $read;
    foreach $read (keys %{$readExtents}) {
	my $extent = $readExtents->{$read};
	my $c1 = $extent->{c1};
	my $c2 = $extent->{c2};
	my $needToTag = &_overlaps($c1, $c2, $leftEnd, $rightEnd);

	if ($needToTag) {
	    my $readLeft = &_posInChild($extent, $leftEnd, 'left');
	    my $readRight = &_posInChild($extent, $rightEnd, 'right');
	    # print STDERR "Need to tag $read $c1, $c2, $leftEnd, $rightEnd $readLeft $readRight\n";
	    if ($readLeft < 1) { $readLeft = 1; }
	    if ($readLeft > $extent->{readLen}) { $readLeft = $extent->{readLen}; }
	    if ($readRight < 1) { $readRight = 1; }
	    if ($readRight > $extent->{readLen}) { $readRight = $extent->{readLen}; }
	    my $areasToTag = new Uncovered($readLeft, $readRight);
	    foreach (split(/\n/, $self->{CAF}->{Sequence}->{$read})) {
		my ($oldtag,$r1, $r2);
		if (($oldtag,$r1, $r2) = /^(?:Tag|Seq_vec|Clone_vec|Clipping)\s+(\S+)\s+(\d+)\s+(\d+)/) {
		    $areasToTag->cover($r1, $r2) if grep($_ eq $oldtag, @{$exclude});
		}
	    }
	    my $area;
	    for($area = $areasToTag->getFirstRegion(); $area; $area = $areasToTag->getNextRegion()) {
		$self->{CAF}->{Sequence}->{$read} .= "$tagKey $newTag $area->[0] $area->[1] \"$comment\"\n";
		# print STDERR "$read $area->[0] $area->[1] $readLeft $readRight $leftEnd $rightEnd\n";
	    }
	}
    }
}

sub _overlaps {
    my ($al, $ar, $bl, $br) = @_;

    if ($al > $ar) { ($al, $ar) = ($ar, $al); }
    if ($bl > $br) { ($bl, $br) = ($br, $bl); }
    
    my $left  = ($al > $bl) ? $al : $bl;
    my $right = ($ar < $br) ? $ar : $br;

    return ($left <= $right);
}

#sub _posInChild {
#    my ($extent, $posInContig) = @_;
#    # my @afLines = sort { $a->{c2} <=> $b->{c2}; } $extent->{afLines};
#    my $afLines = $extent->{afLines};
#    my $last_af = $afLines->[$#$afLines];
#    if ($extent->{c1} < $extent->{c2}) { # Forward strand
#	my $af;
#	foreach $af (@{$afLines}) {
#	    if ($posInContig <= $af->{c2}) {
#		return ($posInContig - $af->{c1} + $af->{r1});
#	    }
#	}
#	return ($posInContig - $last_af->{c1} + $last_af->{r1});
#    } else { # reverse strand
#	my $af;
#	foreach $af (@{$afLines}) {
#	    if ($posInContig <= $af->{c1}) {
#		return ($af->{c2} - $posInContig + $af->{r2});
#	    }
#	}
#	return ($last_af->{c1} - $posInContig + $last_af->{r1});
#    }
#}

sub _posInChild {
    my ($extent, $pos_in_contig, $end) = @_;

    my $af_lines = $extent->{afLines};
    my $last_af = $#$af_lines;

    # Find the af_line for where pos_in_contig is located.
    # Note that for depadded assemblies, pos_in_contig my be at a pad
    # position.  In this case, for the left end of a region, we need the
    # af_line to the right of pos_in_contig, otherwise we need the one to
    # the left.

    my $i;
    if ($end eq 'left') {
	for ($i = 0; $i <= $last_af; $i++) {
	    my $af = $af_lines->[$i];
	    my $rightmost = ($af->{c1} > $af->{c2}) ? $af->{c1} : $af->{c2};
	    if ($rightmost >= $pos_in_contig) { last; }
	}
    } else {
	for ($i = $last_af; $i >= 0; $i--) {
	    my $af = $af_lines->[$i];
	    my $leftmost = ($af->{c1} < $af->{c2}) ? $af->{c1} : $af->{c2};
	    if ($leftmost <= $pos_in_contig) { last; }
	}
    }

    my $af = $af_lines->[$i];

    # Work out where pos_in_contig is in the read.  Note it may be moved
    # slightly so that it doesn't land on a pad.

    if ($extent->{c1} < $extent->{c2}) { # Forward strand
	if ($pos_in_contig < $af->{c1} && $i > 0) {
	    return $af->{r1};  # On a pad to left of current af_line
	} elsif ($pos_in_contig > $af->{c2} && $i < $last_af) {
	    return $af->{r2};  # On a pad to right of current af_line
	} else {
	    # Either in range of current af line, or in cutoff to either side
	    return $pos_in_contig - $af->{c1} + $af->{r1};
	}
    } else { # Reverse strand
	if ($pos_in_contig < $af->{c2} && $i > 0) {
	    return $af->{r2};  # On a pad to left of current af_line
	} elsif ($pos_in_contig > $af->{c1} && $i < $last_af) {
	    return $af->{r1};  # On a pad to right of current af_line
	} else {
	    # Either in range of current af line, or in cutoff to either side
	    return $af->{c2} - $pos_in_contig + $af->{r2};
	}
    }
    die "This statement should be unreachable";
}

sub _getReadExtents {
    my ($self, $contig) = @_;
    croak "Wrong type" unless ref($self);
    my $readExtents = {};
    my $extendLimit = 10;

    foreach (split(/\n/,$self->{CAF}->{Sequence}->{$contig})) {
	my ($read, $c1, $c2, $r1, $r2);
	if (($read, $c1, $c2, $r1, $r2) = /^Assembled_from\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/) {
	    if ($readExtents->{$read}) {
		my $extents = $readExtents->{$read};
		if ($extents->{r1} > $r1) {
		    $extents->{r1} = $r1;
		    $extents->{c1} = $c1;
		}
		if ($extents->{r2} < $r2) {
		    $extents->{r2} = $r2;
		    $extents->{c2} = $c2;
		}
		push(@{$extents->{afLines}}, {r1 => $r1, r2 => $r2, c1 => $c1, c2 => $c2});
	    } else {
		$readExtents->{$read} = {r1 => $r1, r2 => $r2, c1 => $c1, c2 => $c2, readLen => length($self->{CAF}->{DNA}->{$read}),
				        afLines => [{r1 => $r1, r2 => $r2, c1 => $c1, c2 => $c2}]};
	    }
	}
    }
    my $read;
    foreach $read (keys %{$readExtents}) {
	my $extents = $readExtents->{$read};

	# sort assembled_from lines into contig order
	my @afLines = sort { $a->{c2} <=> $b->{c2}; } @{$extents->{afLines}};
	$extents->{afLines} = \@afLines;

#	# extrapolate for cutoffs of readings
#	my $readLen = $extents->{readLen};
#	my $leftAdjust = $extents->{r1} - 1;
#	my $rightAdjust = $readLen - $extents->{r2};
#	$rightAdjust = $extendLimit if ($rightAdjust > $extendLimit);
#	$extents->{r1} = 1;
#	$extents->{r2} += $rightAdjust;
#	if ($extents->{c1} < $extents->{c2}) { # Forward strand
#	    $extents->{c1} -= $leftAdjust;
#	    $extents->{c2} += $rightAdjust;
#	} else { # reverse strand
#	    $extents->{c1} += $leftAdjust;
#	    $extents->{c2} -= $rightAdjust;
#	}
    }
    return $readExtents;
}

sub addTag {
    my ($self, $read, $tagType, $tagStart, $tagEnd, $tagComment) = @_;
    $self->{CAF}->{Sequence} = {} unless $self->{CAF}->{Sequence};
    $self->{CAF}->{Sequence}->{$read} = "" unless defined($self->{CAF}->{Sequence}->{$read});
    if ($tagType eq "SVEC") {
	$self->{CAF}->{Sequence}->{$read} .= "Seq_vec $tagType $tagStart $tagEnd \"$tagComment\"\n";
    } elsif ($tagType eq "CVEC") {
	$self->{CAF}->{Sequence}->{$read} .= "Clone_vec $tagType $tagStart $tagEnd \"$tagComment\"\n";
    } elsif ($tagType eq "STOL") {
	$self->{CAF}->{Sequence}->{$read} .= (defined $tagComment) ? "Stolen \"$tagComment\"\n" : "Stolen\n";
    } else {
	$self->{CAF}->{Sequence}->{$read} .= "Tag $tagType $tagStart $tagEnd \"$tagComment\"\n";
    }
}

sub cleanupTags {
    my ($self, $tagid) = @_;
    confess "Wrong type" unless ref($self);
    my $read;
    foreach $read ($self->readsList(), $self->contigsList()) {
	$self->{CAF}->{Sequence}->{$read} =~ s/(Tag|Seq_vec|Clone_vec|Clipping)\s+$tagid.*\n//g;
    }
}

sub dumpFasta {
    my ($self, $file) = @_;
    confess "Wrong type" unless ref $self;
    my $handle = (ref($file)) ? $file : (FileHandle->new("> $file") || die "Couldn't open $file for writing\n");
    my $contig;
    my $dumped = 0;
    foreach $contig ($self->contigsList()) {
	print $handle ">$contig\n" || die "Error in dumpFasta when witing contig $contig\n";
	my $DNA = $self->{CAF}->{DNA}->{$contig};
	my $len = length($DNA);
	my $i;
	for ($i = 0; $i < $len; $i+=60) {
	    my $line = substr($DNA, $i, 60);
	    $line =~ tr/ACGT\-xX*/acgtn/;
	    print $handle "$line\n" || die "Error in dumpFasta when witing contig $contig\n";
	}
	print $handle "\n" || die "Error in dumpFasta when witing contig $contig\n";
	$dumped++;
    }
    return $dumped;
}

sub readGAP {
    my ($self, $gapdb) = @_;
    confess "Wrong type" unless ref $self;
    my ($project, $version);
    if (($project, $version) = ( $gapdb =~ /(\S+)\.(\S+)/ )) {
	my $GAP2CAF = FileHandle->new("gap2caf -project $project -version $version -preserve -silent |") || die "Could not read GAP database $gapdb\n";
	$self->readCafFile($GAP2CAF);
	$GAP2CAF->close() || die "Error reading GAP database $gapdb: $! $?\n";
    } else {
	die "readGAP: $gapdb is not a GAP database name\n";
    }
}

sub writeGAP {
    my ($self, $gapdb) = @_;
    confess "Wrong type" unless ref $self;
    my ($project, $version);
    if (($project, $version) = ( $gapdb =~ /(\S+)\.(\S+)/ )) {
	my $CAF2GAP = FileHandle->new("| caf2gap -project $project -version $version -preserve -silent") || die "Could not write GAP database $gapdb\n";
	$self->writeCafFile($CAF2GAP);
	$CAF2GAP->close() || die "Error writing GAP database $gapdb: $! $?\n";
    }
    else {
	warn "$gapdb is not a GAP database name\n";
    }
}

sub _warnError {
    my ($msg, $warnOnError, $retval) = @_;
    warn $msg if ($warnOnError);
    return $retval;
}

sub checkForMissingInfo {
    my ($self, $mustHaveBaseQuality, $warnOnError) = @_;
    confess "Wrong type" unless ref $self;
    my $errorCount = 0;
    my $key = 0;
    warn "Checking CAF data for missing entries\n" if ($self->{verbose});
    foreach $key (keys %{$self->{CAF}->{DNA}}) {
	unless ($mustHaveBaseQuality && $self->{CAF}->{BaseQuality}->{$key}) {
	    $errorCount += &_warnError("DNA entry $key has no corresponding base quality entry\n", $warnOnError, 1);
	}
	unless ($self->{CAF}->{Sequence}->{$key}) {
	    $errorCount += &_warnError("DNA entry $key has no corresponding sequence entry\n", $warnOnError, 1);
	}
    }
    foreach $key (keys %{$self->{CAF}->{BaseQuality}}) {
	unless ($self->{CAF}->{DNA}->{$key}) {
	    $errorCount += &_warnError("BaseQuality entry $key has no corresponding DNA entry\n", $warnOnError, 1);
	}
	unless ($self->{CAF}->{Sequence}->{$key}) {
	    $errorCount += &_warnError("BaseQuality entry $key has no corresponding sequence entry\n", $warnOnError, 1);
	}
    }
    foreach $key (keys %{$self->{CAF}->{Sequence}}) {
	unless ($self->{CAF}->{DNA}->{$key}) {
	    $errorCount += &_warnError("Sequence entry $key has no corresponding DNA entry\n", $warnOnError, 1);
	}
	unless ($mustHaveBaseQuality && $self->{CAF}->{BaseQuality}->{$key}) {
	    $errorCount += &_warnError("Sequence entry $key has no corresponding base quality entry\n", $warnOnError, 1);
	}
    }
    return $errorCount;
}

sub checkRawSequences {
    # Checks Sequence entries for raw reads
    my ($self, $warnOnError) = @_;
    confess "Wrong type" unless ref $self;
    my $verbose = $self->{verbose};
    warn "Checking raw sequence entries for missing information\n" if ($verbose);
    my $key;
    my $errorCount;
    my @seqKeys = keys %{$self->{CAF}->{Sequence}};
    my ($count, $total) = (1, scalar(@seqKeys)); 
    foreach $key (@seqKeys) {
	print STDERR "\rChecking entry $count of $total" if ($verbose);
	if ($self->checkRawSequence($key, $warnOnError)) { $errorCount++; }
	$count++;
    }
    print STDERR "\n" if ($verbose);
    return $errorCount;
}

sub checkRawSequence {
    my ($self, $name, $warnOnError) = @_;
    confess "Wrong type" unless ref $self;
    my $entry = $self->{CAF}->{Sequence}->{$name};
    my $retVal = 0;
    unless (ReadTypes->new($name)->getItemInfo('good')) { $retVal += &_warnError("$name does not match the naming convention\n", $warnOnError, 1); }
    unless ($entry =~ /Is_read/) { $retVal += &_warnError("$name is not a read\n", $warnOnError, 2); }
    unless ($entry =~ /SCF_File/) { $retVal += &_warnError("$name does not have an SCF_File entry", $warnOnError, 4); }
    unless ($entry =~ /Template/) { $retVal += &_warnError("$name does not have a template entry\n", $warnOnError, 8); }
    unless ($entry =~ /Insert_size/) { $retVal += &_warnError("$name does not have an insert size entry", $warnOnError, 16); }
    unless ($entry =~ /Primer/) { $retVal += &_warnError("$name does not have a primer entry\n", $warnOnError, 32); }
    unless ($entry =~ /Strand/) { $retVal += &_warnError("$name does not have a strand entry\n", $warnOnError, 64); }
    unless ($entry =~ /Dye/) { $retVal += &_warnError("$name does not have a dye entry\n", $warnOnError, 128); }
    unless ($entry =~ /ProcessStatus/) { $retVal += &_warnError("$name does not have a process status entry\n", $warnOnError, 256); }
    unless ($entry =~ /Asped/) { $retVal += &_warnError("$name does not have an asped date\n", $warnOnError, 512); }
    unless ($entry =~ /Base_caller/) { $retVal += &_warnError("$name does not have a base caller entry\n", $warnOnError, 1024); }
    return $retVal;
}

sub _chomp {
    my ($s) = @_;
    chomp($s);
    return $s;
}

sub insertExp {
    my ($self, $expFile, $read) = @_;
    confess "Wrong type" unless ref $self;
    $read = _chomp($expFile->getEntry("ID")) unless defined($read);
    die "No read name defined\n" unless ($read);
    $self->{CAF}->{Sequence} = {} unless $self->{CAF}->{Sequence};
    $self->{CAF}->{Sequence}->{$read} = "Is_read\n" unless defined($self->{CAF}->{Sequence}->{$read});
    my $seq = \$self->{CAF}->{Sequence}->{$read};
    my $sq = $expFile->getEntry("SQ");
    die "No original sequence data for '$read'\n" unless ($sq);
    $self->{CAF}->{DNA}->{$read} = $sq;
    $self->{CAF}->{DNA}->{$read} =~ tr/-*/n-/;
    my $sqlen = length($sq);
    my $scfFile = &_chomp($expFile->getEntry("LN"));
    $$seq .= "SCF_File $scfFile\n" if $scfFile;
    my $tn = &_chomp($expFile->getEntry("TN"));
    $$seq .= "Template $tn\n" if $tn;
    my $range = &_chomp($expFile->getEntry("SI"));
    if ($range) {
	my ($from, $to);
	$$seq .= "Insert_size $from $to\n" if (($from, $to) = ($range =~ /(\d+)\D+(\d+)/));
    }
    my $ligation = &_chomp($expFile->getEntry("LG"));
    $$seq .= "Ligation_no $ligation\n" if ($ligation);
    my $primer = &_chomp($expFile->getEntry("PR"));
    $$seq .= ($primer == 1 || $primer == 2) ? "Primer Universal_primer\n" : "Primer Custom\n" if ($primer);
    my $strand = &_chomp($expFile->getEntry("DR"));
    $$seq .= ($strand eq "+") ? "Strand Forward\n" : "Strand Reverse\n" if ($strand);
    my $dye = &_chomp($expFile->getEntry("CH"));
    $dye = ("terminator" eq ReadTypes->new($read)->getItemInfo("dye")) unless $dye;
    $$seq .= ($dye == 1) ? "Dye Dye_terminator\n" : "Dye Dye_primer\n";
    my $cloneName = &_chomp($expFile->getEntry("CN"));
    $$seq .= "Clone $cloneName\n" if ($cloneName);
    my $tag;
    for ($tag = $expFile->getFirstEntry("TG"); $tag; $tag = $expFile->getNextEntry("TG")) {
	$tag =~ s/[\n\"]//g;

	# Eg of TG line:
	# TG ALUS + 123..456
	# TG      Alu found by hmmfs

	my ($tagType, $tagDir, $tagStart, $tagEnd, $tagComment) = $tag =~ /^(\S\S\S\S)\s+(\S)\s+(\d+)..(\d+)\s*(.*)$/;
	($tagStart, $tagEnd) = ($tagEnd, $tagStart) if $tagDir eq "-";
	if ($tagType eq "STOL") {
	    $self->addTag($read, $tagType, $tagStart, $tagEnd, $cloneName);
	} else {
	    $self->addTag($read, $tagType, $tagStart, $tagEnd, $tagComment);
	}	
    }
    my $cloneVecPos = &_chomp($expFile->getEntry("CS"));
    if ($cloneVecPos) {
	my ($tagStart , $tagEnd) = $cloneVecPos =~ /(\d+)..(\d+)/;
	my $vector = &_chomp($expFile->getEntry("CV"));
	$$seq .= "Clone_vec CVEC $tagStart $tagEnd $vector\n";
    }
    my $seqVec = &_chomp($expFile->getEntry("SV"));
    my $seqVecLeft = &_chomp($expFile->getEntry("SL"));
    my $seqVecRight = &_chomp($expFile->getEntry("SR"));
    if ($seqVecLeft == $sqlen and $seqVecRight == 1) {
	$$seq .= "Seq_vec SVEC 1 $sqlen $seqVec\n";  # Prevent double entry if completely sequencing vector
    } else {
	$$seq .= "Seq_vec SVEC 1 $seqVecLeft $seqVec\n" if ($seqVecLeft && $seqVecLeft > 0);
	$$seq .= "Seq_vec SVEC $seqVecRight $sqlen $seqVec\n" if ($seqVecRight && $seqVecRight < $sqlen);
    }	
    my $qualLeft = &_chomp($expFile->getEntry("QL"));
    my $qualRight = &_chomp($expFile->getEntry("QR"));
    if ($qualLeft && $qualRight) {
	$qualLeft = 1 if $qualLeft < 1;
	$qualRight = $sqlen if $qualRight > $sqlen;
	$$seq .= "Clipping QUAL $qualLeft $qualRight\n";
    }
    my $processStatus = &_chomp($expFile->getEntry("PS"));
    if ($processStatus) {
	if ($processStatus =~ /Trace quality/i) {
	    $$seq .= "ProcessStatus QUAL  \"$processStatus\"\n";
	} elsif ($processStatus =~ /Matches/i) {
	    $$seq .= "ProcessStatus CONT  \"$processStatus\"\n";
	} elsif ($processStatus =~ /Completely cloning vector/i) {
	    $$seq .= "ProcessStatus CVEC  \"$processStatus\"\n";
	} elsif ($processStatus =~ /Completely sequencing vector/i) {
	    $$seq .= "ProcessStatus SVEC  \"$processStatus\"\n";
	} else {
	    $$seq .= "ProcessStatus FAIL  \"$processStatus\"\n";
	}
    } else {
	$$seq .= "ProcessStatus PASS\n";
    }
    my $aspedDate = &_chomp($expFile->getEntry("DT"));
    $$seq .= "Asped $aspedDate\n" if ($aspedDate);
    my $originalNumbering = &_chomp($expFile->getEntry("ON"));
    if ($originalNumbering) {
	$$seq .= "Padded\n"; 
	my $pos = 1;
	my ($from, $to);
	foreach (split(/\s+/,$originalNumbering)) {
	    if ($_ == 0) {
		$pos++;
	    } elsif (($from,$to) = /^(\d+)\.\.(\d+)/) {
		my $len = $to - $from + 1;
		my $end = $pos + $len - 1;
		$$seq .= "Align_to_SCF $pos $end $from $to\n";
		$pos += $len;
	    } else {
		$$seq .= "Align_to_SCF $pos $pos $_ $_\n";
		$pos++;
	    }
	}
    } else {
	$$seq .= "Unpadded\n";
	$$seq .= "Align_to_SCF 1 $sqlen 1 $sqlen\n";
    }
    my $baseQual = &_chomp($expFile->getEntry("AV"));
    if (defined($baseQual)) {
	my @qual;
	my $val;
	foreach $val (split(/\s+/, $baseQual)) {
	    push(@qual, ($val == 100) ? 99 : $val);
	}
	$self->{CAF}->{BaseQuality}->{$read} = join(" ", @qual);
    }
    my $baseCaller = &_chomp($expFile->getEntry("BC"));
    $$seq .= "Base_caller $baseCaller\n" if ($baseCaller);
    return 0;
}

# Contig version (save consensus tags option in GAP)

sub contigExp2Caf {
    my ($self, $expFile, $read) = @_;
    confess "Wrong type" unless ref $self;
    $self->{CAF}->{Sequence} = {} unless ($self->{CAF}->{Sequence});
    $self->{CAF}->{Sequence}->{$read} = "Is_contig\n" unless (defined($self->{CAF}->{Sequence}->{$read}));
    my $sq = $expFile->getEntry("SQ");
    unless ($sq) {
	warn "No original sequence data for $read\n";
	return;
    }
    my $sqlen = length($sq);
    my $tag;
    for ($tag = $expFile->getFirstEntry("TC"); $tag; $tag = $expFile->getNextEntry("TC")) {
	$tag =~ s/[\n\"]//g;
	my ($tagType, $tagDir, $tagStart, $tagEnd, $tagComment) = $tag =~ /^(\S\S\S\S)\s+(\S)\s+(\d+)..(\d+)\s*(.*)$/;
	($tagStart, $tagEnd) = ($tagEnd, $tagStart) if $tagDir eq "-";
	$self->addTag($read, $tagType, $tagStart, $tagEnd, $tagComment);
    }    
}

sub parse {
    my ($self) = @_;
    confess "Wrong type" unless ref $self;
    my $verbose = $self->{verbose};
    print STDERR "Parsing CAF data\n" if ($verbose);
    my @parsedEntries;
    my @seqKeys = keys %{$self->{CAF}->{Sequence}};
    my ($count, $total) = (0, scalar(@seqKeys));
    my $readName;
    foreach $readName (@seqKeys) {
	$count++;
	print STDERR "\rParsing read $count of $total" if ($verbose);
	my $read = {};
	$read->{readName} = $readName;
	$read->{seqVecs} = [];
	$read->{clippings} = [];
	$read->{tags} = [];
	$read->{cloneVecs} = [];
	$read->{aligns} = [];
	my @lines = split('\n', $self->{CAF}->{Sequence}->{$readName});
	my $line;
	foreach $line (@lines) {
	    if ($line =~ /Template\s*(\S+)/) {
		$read->{template} = $1;
	    } elsif ($line =~ /Strand\s+(Reverse|Forward)/) {
		$read->{strand} = $1;
	    } elsif ($line =~ /Primer\s+(Unknown_primer|Universal_primer|Custom\s*\S*)/) {
		$read->{primer} = $1;
	    } elsif ($line =~ /Dye\s+(Dye_terminator|Dye_primer)/) {
		$read->{dye} = $1;
	    } elsif ($line =~ /Asped\s+(\d\d\d\d-\d\d?-\d\d?)/) {
		$read->{asped} = $1;
	    } elsif ($line =~ /Base_caller\s+(\S+)/) {
		$read->{basecaller} = $1;
	    } elsif ($line =~ /ProcessStatus\s+(\S+)\s*(.*)/) {
		$read->{processStatus} = $1.(($2)?" $2":"");
	    } elsif ($line =~ /SCF_File\s+(\S+)/) {
		$read->{scfFile} = $1;
	    } elsif ($line =~ /Insert_size\s+(\d+)\s+(\d+)/) {
		$read->{insertSize} = { min => $1, max => $2 };
	    } elsif ($line =~ /Ligation_no\s(\S+)/) {
		$read->{ligationNum} = $1;
	    } elsif ($line =~ /Clone\s(\S+)/) {
		$read->{cloneName} = $1;
	    } elsif ($line =~ /Seq_vec\s+SVEC\s+(\d+)\s+(\d+)\s*(\w*)/) {
		push(@{$read->{seqVecs}},
		     { begin => $1, end => $2, text => $3 });
	    } elsif ($line =~ /Clipping\s+(\w+)\s+(\d+)\s+(\d+)\s*(\w*)/) {
		push(@{$read->{clippings}},
		     { type => $1, begin => $2, end => $3, text => $4 });
	    } elsif ($line =~ /Tag\s+(\w+)\s+(\d+)\s+(\d+)\s+\"(.*)\"/) {
		push(@{$read->{tags}},
		     { type => $1, begin => $2, end => $3, text => $4 });
	    } elsif ($line =~ /Tag\s+(\w+)\s+(\d+)\s+(\d+)\s*(\w*)/) {
		push(@{$read->{tags}},
		     { type => $1, begin => $2, end => $3, text => $4 });
	    } elsif ($line =~ /Clone_vec\s+CVEC\s+(\d+)\s+(\d+)\s*(\w*)/) {
		push(@{$read->{cloneVecs}},
		     { begin => $1, end => $2, text => $3 });
	    } elsif ($line =~ /Align_to_SCF\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/) {
		push(@{$read->{aligns}},
		     { startInSeq => $1, endInSeq => $2,
		       startInSCF => $3, endInSCF => $4});
	    }
	}
	$read->{DNA} = $self->{CAF}->{DNA}->{$readName};
	$read->{baseQuality} = pack('C*', split(" ", $self->{CAF}->{BaseQuality}->{$readName}));
	$read->{"length"} = length($read->{DNA});
	$read->{GC} = $read->{DNA} =~ tr/GCgc//;
	push(@parsedEntries, $read);
    }
    print STDERR "\n" if ($verbose);
    return \@parsedEntries;
}

sub addCaf {
    my ($self, $other, $dontOverwrite) = @_;
    confess "Wrong type" unless ref $self;
    my @otherEntries = $other->entriesList();
    my $entry;
    foreach $entry (@otherEntries) {
	next if ($dontOverwrite && defined($self->{CAF}->{Sequence}->{$entry}));
	$self->insertEntry($other->getEntry($entry), $entry);
    }
}

sub process_with {
    my ($self, $tool) = @_;

    $! = 0;
    my $temp = $tFiles->getTemp();
    $self->writeCafFile($temp);
    my $in = FileHandle->new("$tool < $temp |")
	|| die "Couldn't open pipe to $tool $!\n";
    $self->trashAll();
    $self->readCafFile($in);
    close($in) || die "Error $? running $tool $!\n";
    $tFiles->destroyTemp($temp);
}

sub depad {
    my ($self) = @_;
    confess "Wrong type" unless ref $self;
    my $temp = $tFiles->getTemp();
    $self->writeCafFile($temp);
    my $in = FileHandle->new("caf_depad < $temp |") || die "Couldn't open pipe to caf_depad\n";
    $self->trashAll();
    $self->readCafFile($in);
    close($in) || die "Error $? running caf_depad: $!\n";
    $tFiles->destroyTemp($temp);
}

sub pad {
    my ($self) = @_;
    confess "Wrong type" unless ref $self;
    my $temp = $tFiles->getTemp();
    $self->writeCafFile($temp);
    my $in = FileHandle->new("caf_pad < $temp |") || die "Couldn't open pipe to caf_pad\n";
    $self->trashAll();
    $self->readCafFile($in);
    close($in) || die "Error $? running caf_pad: $!\n";
    $tFiles->destroyTemp($temp);
}

sub split {
    my ($self) = @_;
    confess "Wrong type" unless ref $self;
    my $temp = $tFiles->getTemp();
    $self->writeCafFile($temp);
    my $in = FileHandle->new("cafsplit < $temp |") || die "Couldn't open pipe to cafsplit\n";
    $self->trashAll();
    $self->readCafFile($in);
    close($in) || die "Error $? running cafsplit: $!\n";
    $tFiles->destroyTemp($temp);
}

sub autoedit {
    my ($self, $npEditOpts) = @_;
    confess "Wrong type" unless ref $self;
    my $temp = $tFiles->getTemp();
    my $logTemp = $tFiles->getTemp();
    $self->writeCafFile($temp);
    my $in = FileHandle->new("np_edit $npEditOpts -caf $temp -logfile $logTemp | nd_edit | cafcat |") || die "Couldn't open pipe to autoedit\n";
    $self->trashAll();
    $self->readCafFile($in);
    close($in) || die "Error $? running autoedit: $!\n";
    $tFiles->destroyTemp($temp);
    my $log = FileHandle->new("< $logTemp") || die "Couldn't open autoedit log file $logTemp\n";
    my @failures;
    while (<$log>) {
	if ( /misaligned/ ) {
	    my ($read) = split;
	    push(@failures, $read);
	}
    }
    close($log);
    $tFiles->destroyTemp($logTemp);
    return \@failures;
}

sub clip {
    my ($self, $clip, $extend) = @_;
    confess "Wrong type" unless ref $self;
    my @failures;
    my $temp = $tFiles->getTemp();
    $self->writeCafFile($temp);
    $self->trashAll();
    my $svecClip = "nd_clip -caf $temp -svector -logfile /dev/null";
    $clip = 'QUAL' if ($clip eq 'asp');
    if ($clip eq 'eclip') {
	my $in = FileHandle->new("$svecClip | ne_clip -original -logfile /dev/null |") || die "Couldn't open pipe to ne_clip\n";
	$self->readCafFile($in);
	close($in) || die "Error $? running ne_clip: $!\n";
    } else {
	my $logTemp = $tFiles->getTemp();
	my $ext = ($extend) ? "cafextend -clip $clip |" : "";
	my $in = FileHandle->new("$svecClip | $ext nd_clip -clip cafextend,$clip -logfile $logTemp |") || die "Couldn't open pipe to nd_clip\n";
	$self->readCafFile($in);
	close($in) || die "Error $? running nd_clip: $!\n";
	my $log = FileHandle->new("< $logTemp") || die "Couldn't open log file $logTemp\n";
	while (<$log>) {
	    if ( /NDCLIP Failed/ ) {
		my ($read) = split;
		push(@failures, $read);
	    }
	}
	close($log);
	$tFiles->destroyTemp($logTemp);
    }
    $tFiles->destroyTemp($temp);
    return \@failures;
}

sub crossMatch {
    my ($self, $other, $opts) = @_;
    my $ulimit = ($opts->{dlimit}) ? "ulimit -d $opts->{dlimit};" : "";
    my $minscore = ($opts->{minscore}) ? "-minscore $opts->{minscore}" : "";
    my $minmatch = ($opts->{minmatch}) ? "-minscore $opts->{minmatch}" : "";

    my $selfFasta = $tFiles->getTemp();
    print STDERR "Writing caf reads for cross_match (1)\n" if ($self->{verbose});
    my $caf2phrap = FileHandle->new("| cafcat | caf2phrap -fasta $selfFasta -high 1 -low 0 -mask 0 -clipping QUAL -noqual") 
	or die "Couldn't open pipe to caf_phrap\n";
    $self->writeCafFile($caf2phrap);
    close($caf2phrap) or die "Error running caf2phrap: $? $!\n";
    unlink("$selfFasta.caf");

    my $otherFasta = $tFiles->getTemp();
    print STDERR "Writing caf reads for cross_match (2)\n" if ($self->{verbose});
    $caf2phrap = FileHandle->new("| cafcat | caf2phrap -fasta $otherFasta -high 1 -low 0 -mask 0 -clipping QUAL -noqual") 
	or die "Couldn't open pipe to caf_phrap\n";
    $other->writeCafFile($caf2phrap);
    close($caf2phrap) or die "Error running caf2phrap: $? $!\n";
    unlink("$otherFasta.caf");

    print STDERR "Running cross_match...\n" if ($self->{verbose});
    my $cMatch = new FileHandle("$ulimit cross_match $minmatch $minscore $selfFasta $otherFasta |");
    $cMatch or die "Couldn't open pipe to cross_match\n";
    my @cMatchResults;
    while (<$cMatch>) {
	my ($score, $subs, $dels, $inserts, $id1, $start1, $end1, $extra1, $complement, $id2, $start2, $end2, $extra2, $better) = &_readCrossMatch($_);
	if ($score) { push(@cMatchResults, {score => $score, subs => $subs, inserts => $inserts, id1 => $id1, start1 => $start1, end1 => $end1,
					    extra1 => $extra1, id2 => $id2, start2 => $start2, end2 => $end2, extra2 => $extra2, better => $better}); }
    }
    close($cMatch) or die "cross_match returned error $?: $!\n";
    unlink("$selfFasta.log");
    $tFiles->destroyTemp($selfFasta, $otherFasta);
    return \@cMatchResults;
}

sub _readCrossMatch {
    my ($line) = @_;
    chomp($line);
    my ($score, $subs, $dels, $inserts, $id1, $start1, $end1, $extra1, $complement, $id2, $extra2, $start2, $end2, $better);
    $line =~ tr/\(\)//d;
    if ($line =~ /\sC\s/) {
	($score, $subs, $dels, $inserts, $id1, $start1, $end1, $extra1, $id2, $extra2, $start2, $end2, $better) = $line =~ /^\s*(\d+)\s+(\d+\.?\d*)\s+(\d+\.?\d*)\s+(\d+\.?\d*)\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)\s+C\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)\s*(\*?)\s*$/;
	$complement = 1;
    } else {
	($score, $subs, $dels, $inserts, $id1, $start1, $end1, $extra1, $id2, $start2, $end2, $extra2, $better) = $line =~ /\s*(\d+)\s+(\d+\.?\d*)\s+(\d+\.?\d*)\s+(\d+\.?\d*)\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)\s*(\*?)\s*$/;
	$complement = 0;
    }
    return ($score, $subs, $dels, $inserts, $id1, $start1, $end1, $extra1, $complement, $id2, $start2, $end2, $extra2, $better);
}

sub assemblyDelReads {
    # Note - assembly may have to go through cafsplit after doing this
    my ($self, @readsToDel) = @_;
    my $verbose = $self->{verbose};
    my %toDel = map { $_, 1 } @readsToDel;
    my $read;
    foreach $read (keys %toDel) {
	$self->trashRead($read);
    }
    my @contigsToDel;
    my $contig;
    foreach $contig ($self->contigsList()) {
	my $contigSeq = \$self->{CAF}->{Sequence}->{$contig};
	my @lines = split(/\n/, $$contigSeq);
	my $assembledFroms = 0;
	$$contigSeq = "";
	foreach (@lines) {
	    if (/Assembled_from\s+(\S+)/) {
		if (!exists($toDel{$1})) {
		    $$contigSeq .= "$_\n";
		    $assembledFroms++;
		}
	    } elsif (/Golden_path\s+(\S+)/) {
		if (!exists($toDel{$1})) { $$contigSeq .= "$_\n"; }
	    } else {
		$$contigSeq .= "$_\n";
	    }
	}
	unless ($assembledFroms) {
	    if ($verbose) {
		print STDERR "All reads from contig $contig have been deleted\n";
	    }
	    push(@contigsToDel, $contig);
	}
    }
    foreach (@contigsToDel) { $self->trashRead($_); }
}

END {
    $tFiles->destroyTemp() if ($tFiles);
}

1;

# $Log: Caftools.pm,v $
# Revision 1.13  2001/06/29  11:12:56  mng
# Added trashContig.
#
# Revision 1.12  2001/04/24  14:12:18  rmd
# Bug fix to tag routine.  Did not work correctly when tagging depadded
# sequences.
#
# Revision 1.11  2001/03/08  11:16:36  rmd
# Changed basequality split in parse function to ignore leading whitespace
#
# Revision 1.10  2001/01/31  09:53:23  rmd
# Added process_with subroutine
#
# Revision 1.9  2000/10/20  13:21:17  mng
# quoted {"length"} to remove the warning message.
#
# Revision 1.8  1999/11/16  11:48:25  rmd
# Added get_entry_ref plus some mods to raw caf parsing function
#
# Revision 1.7  1998/06/12 11:42:25  rmd
# Added method assemblyDelReads - deletes reads from an assembly, and removes
# all references to them from contigs.
#
# Revision 1.6  1998/02/03  16:56:56  sd
# Use ReadTypes instead of a pattern.
#
# Revision 1.5  1997/11/20  14:25:51  rmd
# General improvements
# Fixed bugs in _posInChild
# Support in tag for caching read extents to speed it up
# Added methods to call caf_pad, caf_depad etc.
#
# Revision 1.4  1997/06/20  13:25:02  rmd
# Fixed double Seq_vec entry when completely sequencing vector
#
# Revision 1.3  1997/06/20  10:35:44  rmd
# Fixed bug with readCafFile that prevented it from reading a pipe.
#
# Revision 1.2  1997/06/16  10:35:14  rmd
# Fixed regular expression for process status in parse.
#
# Revision 1.1  1997/06/09  15:43:47  rmd
# Initial revision
#
#
