#! /usr/local/bin/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.                      #
#######################################################################
#
# 13-Feb-1995 Created -- RD
# 14-Feb-1995 Added output to file to writeAceFile -- SD
# 14-Feb-1995 Added &revCompSeq() -- SD
# 15-Feb-1995 Supress empty opjects in writeAceFile -- RD
# 15-Feb-1995 More efficient revCompSeq() -- RD
# 17-Feb-1995 Append lines to class object -- RD
# 22-Feb-1995 eof() check in readAceFile was on STDIN not $file -- SD
# 09-Mar-1995 local($x);
# 15-Mar-1995 Added subroutines theContigs and theReads -- SD
# 10-Apr-1995 Added subroutine trashRead -- SD
# 27-Apr-1995 Tried to improve efficiency of writeAceFile -- SD
# 29-Feb-1996 Support BaseQuality -- SD
# 02-Apr-1996 Renamed to caftools.pl, so as not to conflict with old acetools -rmott
# 03-Apr-1996 Base quality now stored as strings rather than arrays for effienct -- DMH 
#
#
# $Id: caftools.pl,v 1.33 2006/08/29 09:11:48 jkb Exp $
# 
# $Log: caftools.pl,v $
# Revision 1.33  2006/08/29  09:11:48  jkb
# Added use of cover.pl (moved code out of old exptools.pl)
#
# Revision 1.32  1997/04/29 16:56:14  sd
# writeAceFile now dies if error when writing.
#
# Revision 1.31  1997/04/29  16:25:00  sd
# Added a warning to dumpfasta if error during writing.
#
# Revision 1.30  1997/04/24  15:02:02  dn1
# Added copyright
#
# Revision 1.29  1997/02/21  15:13:58  sd
# Added trashAll;
#
# Revision 1.28  1996/11/28  15:01:39  sd
# Initialise BaseQuality to avoid perl -w warnings.
#
# Revision 1.27  1996/11/19  10:41:14  badger
# parenthesised my list.
#
# Revision 1.26  1996/10/17  13:12:13  rmott
# dittot
#
# Revision 1.25  1996/10/17  13:10:29  rmott
# diito
#
# Revision 1.24  1996/10/17  13:09:42  rmott
# ditto
#
# Revision 1.23  1996/10/17  13:04:59  rmott
# ditto
#
# Revision 1.22  1996/10/17  13:01:40  rmott
# ditto
#
# Revision 1.21  1996/10/17  13:00:32  rmott
# added forward-reverse info to readReport
#
# Revision 1.20  1996/10/07  11:39:10  rmott
# ditto
#
# Revision 1.19  1996/10/07  11:36:51  rmott
# changes to read_report
#
# Revision 1.18  1996/10/07  10:37:36  rmott
# added parsing of ProcessStatus to readReport
#
# Revision 1.17  1996/09/30  17:55:09  rmott
# Added declaration for $cloning_vector.
#
# Revision 1.16  1996/09/30  17:51:58  rmott
# A few bugs fixed.
#
# Revision 1.15  1996/09/30  17:46:22  rmott
# ditto ditto
#
# Revision 1.14  1996/09/30  17:44:12  rmott
# ditto
#
# Revision 1.13  1996/09/30  17:41:06  rmott
# added function readReport
#
# Revision 1.12  1996/09/18  10:50:31  rmott
# changed calling syntax for caf2gap & gap2caf
#
# Revision 1.11  1996/09/18  10:32:52  rmott
# bug fix in readGAP
#
# Revision 1.10  1996/09/17  12:58:23  rmott
# added functions readGAP() and writeGAP() for transparent GAP interface
#
# Revision 1.9  1996/07/26  17:04:19  sd
# A number of patterns in parseAceLine were not rooted at the
# start of the line.
#
# Revision 1.8  1996/06/04  09:48:14  sd
# Added parseAceLine(), which splits line into fields, taking
# into account double quotes, etc.
#
# Revision 1.7  1996/05/24  13:37:06  sd
# Added theAssemblies()
#
# Revision 1.6  1996/05/23  08:45:17  sd
# Preserve strand for tags when tagging on readings.
#
# Revision 1.5  1996/05/01  09:48:27  badger
# Sorted comment leader once and for all!
#
# Revision 1.4  1996/05/01  09:43:55  badger
# Same problem with commenting in RCS.
#
# Revision 1.3  1996/05/01  09:42:37  badger
# invalid comment
#
# Revision 1.2  1996/04/10  16:05:13  badger
# handles translation of non-nucleotides correctly - rmott
#
# Revision 1.1  1996/04/03  18:18:01  dmh
# Initial revision
#
#
#
###################################################

require "cover.pl";

sub readAceFile {
#
# read a .ace file into %Class - remove whitespace from DNA
#
    local($file) = @_ ;
    local(@x);
    while (<$file>) {
        s/\/\/.*// ;			# ace file comments
        if ( /^DNA\s+(:\s+)?(\S+)/ ) {
            $name = $2 ;
	    @x = ();
	    while (<$file>) {
		last if /^\s*$/ ;
		s/\s+//g ;		# remove whitespace
		push(@x,$_);
	    }
	    $DNA{$name} = join("",@x) ;
	} elsif ( /BaseQuality\s+(:\s+)?(\S+)/ ) {
            $name = $2 ;
	    $BaseQuality{$name} = "";
	    while (<$file>) {
		last if /^\s*$/ ;
		$BaseQuality{$name} .= $_;
	    }
	} elsif ( /^(\S+)\s+(:\s+)?(\S+)/ ) {
	    $class = $1 ;
	    $name = $3 ;
	    @x = ();
	    while (<$file>) {
		last if /^\s*$/ ;
		s/\\\n$// ;		# continuation lines
		push(@x,$_);
	    }
	    if (defined($$class{$name})) {
	        $$class{$name} .= join("",@x);
	    } else {
	        $$class{$name} = join("",@x);
	    }
	}
	if (eof($file)) { last ; }
    }
}

###################################################

sub writeAceFile {
#
# write a .ace file to first argument
#
    local($file) = @_ ;
    local($key);

    foreach $key (keys %DNA) {
	next unless $DNA{$key};
	die "$prog: Error in writeAceFile while writing CAF file" unless print $file "\nDNA : ",$key,"\n";
        $len = length ($DNA{$key});
	$DNA{$key} =~ s/[^ACGTUacgtu\-\*]/N/g; # translate all non-nucleotides to N's
#	$DNA{$key} =~ s/\*/\-/g; # translate old-style pads
	for ($i = 0 ; $i < $len ; $i += 60) {
	    die "$prog: Error in writeAceFile while writing CAF file" unless print $file substr ($DNA{$key}, $i, 60), "\n";
	}
    }

    foreach $key (keys %BaseQuality) {
	next unless $BaseQuality{$key};
	die "$prog: Error in writeAceFile while writing CAF file" unless print $file "\nBaseQuality : ",$key,"\n";
	die "$prog: Error in writeAceFile while writing CAF file" unless print $file "$BaseQuality{$key}\n";
    }

    foreach $key (keys %Sequence) {
	next unless $Sequence{$key};
	die "$prog: Error in writeAceFile while writing CAF file" unless print $file "\nSequence : ",$key,"\n";
	die "$prog: Error in writeAceFile while writing CAF file" unless print $file $Sequence{$key};
    }


}

###################################################

sub parseAceLine {
    local($line) = @_;
    local(@fields) = ();
    $line =~ /^\s*/;     # Consume opening space
    while ($line = $') {
        if ($line =~ /^"((([^"\\])|(\\.))*)"\s*/ || $line =~ /^(\S+)\s*/ || $line =~ /^.*/) { #"
            push(@fields,$1);
        }
    }
    @fields;
}

###################################################

sub revCompSeq { # -- SD
# 
# Reverse complement a sequence
#
    local($seq) = @_ ;
    $seq = reverse ($seq) ;
    $seq =~ tr/acgtACGT/tgcaTGCA/ ;
    $seq ;
}

###################################################

sub theContigs {
    #
    # Returns a list of contigs
    #
    grep($Sequence{$_} =~ /Is_contig\n/,keys %Sequence);
}

###################################################

sub theReads {
    #
    # Returns a list of readings
    #
    grep($Sequence{$_} =~ /Is_read\n/,keys %Sequence);
}

###################################################

sub theAssemblies {
    #
    # Returns a list of assembly Sequence objects
    #
    grep($Sequence{$_} =~ /Is_assembly\n/,keys %Sequence);
}


###################################################

sub trashRead {
    #
    # Delete a reading from relevant structures
    #
    local( $read ) = @_;
    undef $Sequence{$read};
    undef $DNA{$read};
    undef $BaseQuality{$read};

}


sub trashAll {
    #
    # Delete a reading from relevant structures
    #
    undef %Sequence;
    undef %DNA;
    undef %BaseQuality;

}


###################################################



############################################################
# &tag($contig,$left,$right,$newtag,$comment,\@exclude);
# 
# NOTES:
# [1] Tags readings in region with tag type $newtag
# [2] Do not tag readings where they are already tagged with any in @exclude
############################################################

sub tag {

    local ($contig, $left, $right,$newtag, $comment, $exclude) = @_;
    local (@Af);
    local ($read, $c1, $c2, $r1, $r2);
    local ($i,$j,$extendLimit,$adjust);
    local ($readlen);
    local (@extent);
    local ($key);
    local ($leftend, $rightend);
    
    warn "Need to tag reads for $newtag in $contig from $left through $right\n" if $debug;

    if ($left <= $right) {
	($leftend, $rightend) = ($left, $right);
    } else {
	($leftend, $rightend) = ($right, $left);
    }

    #############################################################
    # Tag consensus
    #############################################################
    if ($newtag eq "CVEC") {
	$Sequence{$contig} .= "Clone_vec $newtag $left $right \"$comment\"\n";
    } else {
	$Sequence{$contig} .= "Tag $newtag $left $right \"$comment\"\n";
    }


    #############################################################
    # Tag readings - a long and involved process
    #############################################################


    #############################################################
    # Unpack contig object
    #############################################################
    @Af = ();

    foreach ( split(/\n/,$Sequence{$contig}) ) {
	
	if (($read, $c1, $c2, $r1, $r2) = /^Assembled_from\s+(\S+)\s+(\d+)\s+(\d+)\s+(\d+)\s+(\d+)/) {
	    push(@Af,[$read,$c1,$c2,$r1,$r2]);
	}
    }

    ############################################################
    # Sort - this may not be necessary but do it anyway
    ############################################################
    @Af = sort AfCmp @Af;


    #############################################################
    # Determine region covered by each reading in contig
    #############################################################

    $extendLimit = 10;

    for ($i = 0; $i <= $#Af; ) {

	# $i points to start of first Af line for the reading

	
	# Determine the extent of coverage of each reading
	@extent = @{$Af[$i]};
	$j = $i;
	for ($i++; $i <= $#Af && $Af[$i]->[0] eq $Af[$j]->[0]; $i++) {
	    if ($extent[3] > $Af[$i]->[3]) {
		$extent[3] = $Af[$i]->[3];
		$extent[1] = $Af[$i]->[1];
	    }
	    if ($extent[4] < $Af[$i]->[4]) {
		$extent[4] = $Af[$i]->[4];
		$extent[2] = $Af[$i]->[2];
	    }
	}

	# Extrapolate for cutoffs of readings
	$readlen = length($DNA{$extent[0]});
	if ($extent[1] < $extent[2]) {
	    # Forward strand
	    # Left end of reading
	    $adjust = $extent[3] - 1;
	    $extent[1] -= $adjust;
	    $extent[3] -= $adjust;
	    # Right end of reading
	    $adjust = $readlen - $extent[4];
	    if ($adjust > $extendLimit) { $adjust = $extendLimit; }
	    $extent[2] += $adjust;
	    $extent[4] += $adjust;
	    
	} else {
	    # Reverse strand
	    # Left end of reading
	    $adjust = $extent[3] - 1;
	    $extent[1] += $adjust;
	    $extent[3] -= $adjust;
	    # Right end of reading
	    $adjust = $readlen - $extent[4];
	    if ($adjust > $extendLimit) { $adjust = $extendLimit; }
	    $extent[2] -= $adjust;
	    $extent[4] += $adjust;
	}

	# If reading overlaps range then we need to tag.
	if ($extent[1] < $extent[2]) {
	    $needToTag = ($extent[2] >= $leftend && $extent[1] <= $rightend);
	} else {
	    $needToTag = ($extent[1] >= $leftend && $extent[2] <= $rightend);
	}


	if ($needToTag) {
	    local(@uncovered);
	    local($readLeft,$readRight);

	    # determine bases in reading requiring tag
	    $readLeft = &posInChild(\@Af,$extent[0],$left);
	    $readRight = &posInChild(\@Af,$extent[0],$right);

	    # clip to reading coords
	    if ($readLeft < 1) { $readLeft = 1; }
	    if ($readLeft > $readlen) { $readLeft = $readlen; }
	    if ($readRight < 1) { $readRight = 1; }
	    if ($readRight > $readlen) { $readRight = $readlen; }

	    # ignore regions tagged with any of @exclude tags (mutually exclusive)
	    &initcover($readLeft,$readRight);
	    foreach ( split(/\n/,$Sequence{$extent[0]}) ) {
		if (($key,$oldtag,$r1, $r2) = /^(Tag|Seq_vec|Clone_vec|Clipping)\s+(\S+)\s+(\d+)\s+(\d+)/) {
		    if (grep($_ eq $oldtag, @{$exclude})) {
			&cover($r1,$r2);
		    }
		}
	    }

	    # Create new Tag entries
	    for ($j=0;$j <= $#uncovered; $j++) {
		if ($newtag eq "CVEC") {
		    if ($readLeft < $readRight) {
			$Sequence{$extent[0]} .= "Clone_vec $newtag $uncovered[$j]->[0] $uncovered[$j]->[1] \"$comment\"\n";
		    } else {
			$Sequence{$extent[0]} .= "Clone_vec $newtag $uncovered[$j]->[1] $uncovered[$j]->[0] \"$comment\"\n";
		    }
		} else {
		    if ($readLeft < $readRight) {
			$Sequence{$extent[0]} .= "Tag $newtag $uncovered[$j]->[0] $uncovered[$j]->[1] \"$comment\"\n";
		    } else {
			$Sequence{$extent[0]} .= "Tag $newtag $uncovered[$j]->[1] $uncovered[$j]->[0] \"$comment\"\n";
		    }
		}
	    }
	}
    }

}




############################################################
# $origpos = &posInChild(\@Af,$read,$pos);
#
# Determines the coordinate of a location in a Sequence
# in the child reading's coordinates.
#
############################################################
sub posInChild {
    local ($Af,$child, $pos) = @_;
    local ($i);
    local ($from,$to);

    # identify subset of Af pertaining to child
    $from = -1;
    $to = -1;
    for ($from = 0; $from <= $#{$Af} && $Af->[$from]->[0] ne $child; $from++) {
    }
    if ($from > $#{$Af}) {
	#not found
	return $pos;
    }
    for ($to = $from + 1; $to <= $#{$Af} && $Af->[$to]->[0] eq $child; $to++) {
    }
    $to--;

    ############################################################
    # Determine orientation
    ############################################################
    if ($Af->[$from]->[1] < $Af->[$to]->[2]) {

	############################################################
	# Find position - FORWARD STRAND
	############################################################

	for($i=$from;$i <= $to; $i++) {
	    if ( $pos <= $Af->[$i]->[2] ) {
		return $pos - $Af->[$i]->[1] + $Af->[$i]->[3];
	    }
	}

	# Still unmapped
	return $pos - $Af->[$to]->[2] + $Af->[$to]->[4];
    } else {

	############################################################
	# Find position - REVERSE STRAND
	############################################################

	for($i=$to;$i >= $from; $i--) {
	    if ( $pos <= $Af->[$i]->[1] ) {
		# OK
		return $Af->[$i]->[2] - $pos + $Af->[$i]->[4];
	    }
	}

	# Still unmapped
	return $Af->[$from]->[1] - $pos + $Af->[$from]->[3];
    }

}


############################################################
# &cleanupTags()
# Removes tags from Readings and Contigs
############################################################

sub cleanupTags {
    local($tagid) = @_;
    local($read);

    # Clean up contigs
    foreach $read ( &theReads , &theContigs ) {
	$Sequence{$read} =~ s/(Tag|Seq_vec|Clone_vec|Clipping)\s+$tagid.*\n//g;
    }

}


############################################################
# $Ncontigs = &dumpfasta(HANDLE)
# Dump consensus sequences for contigs read in using &readAceFile()
# to a handle in fasta format. Returns number of contigs dumped.
############################################################

sub dumpfasta {
    local ($handle) = @_;
    local ($dumped) = 0;
    local ($contig,$line,$i);
    local ($ok);

    foreach $contig ( &theContigs ) {
	print $handle ">$contig\n";
	for($i=0;$i<length($DNA{$contig});$i+=60) {
	    $line = substr($DNA{$contig},$i,60);
	    # CAF to fasta sequence
	    #$line =~ tr/ACGTxX*/acgtn/;
	    $line =~ tr/ACGT\-xX*/acgtn/;
	    $ok = print $handle $line,"\n";
	    die "$0: Error in dumpfasta while writing contig $contig\n" if (! $ok);
	}
	$dumped++;
    }
    $dumped;
}

############################################################
# $length = &seqcount(HANDLE)
# Returns the length of the fasta file read in on handle
############################################################

sub seqcount {
    local ($handle) = @_;
    local($count) = 0;

    while (<$handle>) {
	next if /^>/;
	s/\s//g;
	$count += length($_);
    }

    $count;

}

############################################################
# ($name, $sequence) = &readfasta(HANDLE)
# Reads in a fasta sequence from a handle.
# Name is the first word on the fasta description line.
############################################################

sub readfasta {
    local ($handle) = @_;
    local($seq) = "";
    local($name) = "";

    while (<$handle>) {
	if (/^>\s*(\S*)/) {
	    $name = $1;
	    next;
	}
#	next if (($name) = /^>\s*(\S*)/);
	s/\s//g;
	$seq .= $_;
    }

    ($name,$seq);

}

#############################################################
# &AfCmp()
# Specialised sort subroutine for sorting stacks of Assembled_from arrays.
#############################################################

sub AfCmp {
    local($compare);

    $compare = $a->[0] cmp $b->[0];
    if ($compare==0) {
	if ( ($a->[1] <=> $b->[1]) == ($a->[3] <=> $b->[3]) ) {
	    # forward strand
	    $compare = $a->[1] <=> $b->[1];
	} else {
	    # reverse strand... gotta reverse
	    $compare = $b->[1] <=> $a->[1];
	}
    }

    return $compare;
}






1 ;	# required return value

#############################################################
# &readGAP()
# Read a GAP database via gap2caf
# calling prcedure: &readGAP($gapdb)
# where $gapdb is a string like T06C12.0 giving the name of the database
# Dies if gap2caf won't run
# data are slurped into the standard global structures
#############################################################

sub readGAP {
    my $gapdb = $_[0];
    my $project = "";
    my $version = "";
    if ( ($project,$version) = ( $gapdb =~ /(\S+)\.(\S+)/ ) ) {
	open ( GAP2CAF , "gap2caf -project $project -version $version -preserve |" ) || die "Could not read GAP database $gapdb\n";
	readAceFile( GAP2CAF );
	close( GAP2CAF );
    }
    else {
	warn "$gapdb is not a GAP database name\n";
    }
    0;
}

#############################################################
# &writeGAP()
# writes the current CAF into caf2gap to create a GAP database
# call as &writeGAP($gapdb)
# $gapdb MUST split into PROJECT.VERSION
#############################################################

sub writeGAP {
    my $gapdb = $_[0];
    my $project = "";
    my $version = "";
    if ( ($project,$version) = ( $gapdb =~ /(\S+)\.(\S+)/ ) ) {
	open ( CAF2GAP, "| caf2gap -project $project -version $version -preserve" ) || die "Could not write GAP database $gapdb\n";
	writeAceFile( CAF2GAP );
	close( CAF2GAP );
    }
    else {
	warn "$gapdb is not a GAP database name\n";
    }
    0;
}

#############################################################
# readReport(FILE)
# makes a report on the reads to FILE
#############################################################

sub readReport {
    my ($file) = $_[0];
    my $reads = 0;
    my $full_length = 0;
    my $cloning_vector = 0;
    my $total_len = 0;
    my $total_slen = 0;
    my $reverse = 0;
    my $forward = 0;
    my $fr_pairs = 0;
    my $template;
    my %is_template;
    my %forwards;
    my %reverses;
    my $primers = 0;
    my $terminators = 0;

    foreach $read ( &theReads ) {
	my ($len, $slen, $sl, $sr, $ql, $qr, $cs);
	$reads++;
	$len = length($DNA{$read});
	$sl = 0; $sr = $len + 1;
	$ql = 0; $qr = $len + 1;
	$cs = 0;

	undef($ps);undef($reason);
	undef($template); undef($strand);

	foreach ( split( /\n/, $Sequence{$read} ) ) {
	    if (/^Seq_vec\s+\S+\s+(\d+)\s+(\d+)/) {
		if ($1 == 1) {
		    $sl = $2;
		} elsif ( $2 == $len ) {
		    $sr = $1;
		}
	    }
	    elsif (/^Clone_vec/) {
		$cs++;
	    }
	    elsif (/^Clipping\s+QUAL\s+(\d+)\s+(\d+)/) {
		($ql,$qr) = ($1,$2);
	    }
	    elsif ( /^ProcessStatus\s+(\S+)/ ) {
		$ps = $1;
	    }
	    elsif (/^Strand\s+(\S+)/ ) {
		$strand = $1;
	    }
	    elsif( /^Template\s(\S+)/ ) {
		$template = $1;
	    }
	    elsif( /Dye_terminator/ ) {
		$terminators++;
	    }
	    elsif (/Dye_primer/ ) {
		$primers++;
	    }
	}

	$full_length++ if ($sr != $len + 1);
	$cloning_vector ++ if $cs;
	$total_len += $len;
	$total_slen += $sr - $sl - 1;
	$total_clen += $qr - $ql + 1;
	$PS{$ps}++;
	$CONT{$reason}++ if ( $ps eq "CONT" );
	if ( $strand eq "Forward" ) {
	    $forward ++;
	    $forwards{$template}++;
	}
	if ( $strand eq "Reverse" ) {
	    $reverse ++ ;
	    $reverses{$template}++;
	}
	push @{$reads_off{$template}}, $strand;
	$is_template{$template}++;
    }

    $templates = scalar (keys %is_template);
    foreach $template (keys %is_template) {
	$fr_pairs++ if ( $forwards{$template} && $reverses{$template} );
    }

    print $file "Number of reads                           = $reads\n";
    print $file "Number of primer reads                    = $primers\n";
    print $file "Number of terminator reads                = $terminators\n";
    print $file "Number of forward reads                   = $forward\n";
    print $file "Number of reverse reads                   = $reverse\n";
    print $file "Number of templates                       = $templates\n";
    print $file "Number of forward-reverse pairs           = $fr_pairs\n";
    print $file "Number of full length reads               = $full_length\n";
    print $file "Number of reads containing cloning vector = $cloning_vector\n";

    printf $file "Average read length                       = %6.1f\n", $total_len / $reads;
    printf $file "Average read length (vector clipped)      = %6.1f\n", $total_slen / $reads;
    printf $file "Average read length (quality clipped)     = %6.1f\n", $total_clen / $reads;
    print "Process Status:\n";
    foreach $ps (keys %PS) {
	printf "%-10s %5d\n", $ps, $PS{$ps};
    }
}

############# end of file ###############
# Local Variables:
# mode:perl
# End:
############# end of file ###############
