#! /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.                      #
#######################################################################
#
# 22-Feb-1995 Created -- SD
# 08-Mar-1995 Added some support for staden format files -- SD
# 09-Mar-1995 Bug resulted in last line of experiment file being dropped -- SD
# 05-Apr-1995 Added &cover() -- SD
# 19-Apr-1995 readExpFile() now returns status -- SD
# 28-Apr-1995 Added &killEntry and &replaceEntry() -- DMH
# 27-Jun-1995 Cleaned up some -w switch warnings -- SD
# 27-Jun-1995 Bug in &initcover() if $c1 > $c2 -- SD
# 09-Aug-1995 Serious bug! Was using %Nentries instead of %expNentries in a few places -- SD
# 06-Feb-1996 Bug enforcing terminal newline in &addEntry() -- SD
# 15-Feb-1996 Added special treatment for writing AV field -- MCW
# 01-Apr-1996 Conversion of - to N in DNA -- RMOTT
#


############################################################
# err = &readExpFile ( FILEHANDLE )
#
# Sets global variables
#    @expOrder
#	 Order in which exp entries occur
#    %expNentries
#        Number of lines for each type of entry
#        eg $expNentries{"TG"} is the number of tag lines
#    %expEntry
#        Array of lines for each entry
#        eg $expEntry{"TG"}[0] is the first tag line
#        eg $expEntry{"TG"}[$expNentries{"TG"}] is the last tag line
#
# NOTES on entries:
# [1] Lines contain \n on the end. May need to chop
# [2] Multiline entries (eg TG) appear in one entry, lines separated by  \n
# [3] Sequence (SQ) has white space removed for convenience
# [4] RETURNS error code. 0 means no error.
#
############################################################

sub readExpFile {
    local($file) = @_ ;
    local($multiline);
    local($id);
    local($n);
    local($value);

    undef @expOrder;
    undef %expNentries;
    undef %expEntry;

    # read the first line
    $_ = <$file>;

    if (($id) = /^ID\s+(\S+)/) {
	############################################################
	# experiment file
	############################################################
	addEntry("ID",$id);
	$multiline = "";
	while (<$file>) {
	    ( $id ) = /^(..)/;
	    if ($id =~ /^\s+$/ && $multiline ne "") {
		$id = $multiline;
	    }
	    if ($id eq "SQ") {
		############################################################
		# special treatment for sequences
		############################################################
		push(@expOrder,$id); #keep a track of order
		$n = $expNentries{$id};
		$multiline = "";
		while (<$file>) {
		    last if /^\/\//;
		    $_ =~ s/\s//g; # Kill white space
#		    $_ =~ s/\-/N/g; # Replace '-' by 'N' 
		    $expEntry{$id}[$n] .= $_;
		}
		$expNentries{$id}++;
	    } else {
		$value = substr($_,5);
		$n = $expNentries{$id}+0;
		if ($id eq $multiline && $value =~ /^     /) {
		    $expEntry{$id}[$n-1] .= $value;
		} else {
		    push(@expOrder,$id); #Keep a track of line order
		    $multiline = "";
		    $expEntry{$id}[$n] = $value;
		    $expNentries{$id}++;
		    if ($id eq "TG" || $id eq "TC" || $id eq "ON" || $id eq "AV") { $multiline = $id; }
		}
	    }
	}
    } elsif (/^;/) {
	############################################################
	# staden format file
	############################################################
	local($left,$right,$seq,$type,$from,$to,$len,$comment);
	local($leftlen,$rightlen,$seqlen);
	($lt,$ln) = /;\s*\d+\s*\d+\s*\d+(....)(\S+)\s*$/;
	$lt =~ s/ //g;
	addEntry("LT",$lt);
	addEntry("LN",$ln);
	while (<$file>) {
	    if (/^;</) {
		$left .= substr($_,2);
	    } elsif (/^;>/) {
		$right .= substr($_,2);
	    } elsif (/^;;/) {
		#;;ALUS      1    146      0.873
		($type,$from,$len,$comment) = /^;;(....)\s+(\d+)\s+(\d+)\s*(.*)$/;
		$to = $from + $len - 1;

		# FIXME YUK YIKES BUG HACK
		# BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG
	        # In staden format, tag positions are relative to the first used base (less one)
	        # and must not fall outside the reading.
                # Currently this is not supported
	        # BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG BUG

		if (defined($comment)) {
		    addEntry("TG","$type + $from..$to\n     $comment");
		} else {
		    addEntry("TG","$type + $from..$to");
		}
	    } elsif (/^;/) {
		# Comment
		addEntry("CC",substr($_,1));
	    } else {
		$seq .= $_;
	    }
	}

	# remove white space
	$left =~ s/\s//g;  $leftlen = length($left);
	$right =~ s/\s//g; $rightlen = length($right);
	$seq =~ s/\s//g;   $seqlen = length($seq);
	#Now dump out sequence...
	$seq = $left . $seq . $right;
	addEntry("QL",$leftlen);
	addEntry("QR",$leftlen+$seqlen+1);
	addEntry("SQ",$seq);
    } elsif (/^>(\S)/)  { # fasta sequence file
	addEntry("ID",$1);
	undef $seq;
	while(<$file>) {
	    chop();
	    s/\s*//;
	    $seq .= $_;
	}
	if ( defined($seq) ) {
	    addEntry("SQ", $seq );
        }
    } elsif (/^[\s\*ACGT\-acgtNn]+$/)  { # plain sequence file
	s/\s*//g;
	$seq = $_;
	while(<$file>) {
	    chop();
	    if  (/^[\s\*ACGT\-acgtNn]+$/) {
		s/\s*//g;
		$seq .= $_;
	    } else {
		undef($seq);
		last;
	    }
	}
	if ( defined($seq) ) {
	    addEntry("SQ", $seq );
        } else {
	    warn "readExpFile: unknown format - cannot read $file\n";
	    return 1;
        }
    } else {		
        warn "readExpFile: unknown format - cannot read $file\n";
        return 1; 
    }
    0;				
}

############################################################
# &getEntry ( $lineType )
#
# Returns the last line of specified line type
# eg &getEntry("QL") returns the last quality clipping line
# if there is no entry, it returns undefined.
# 
############################################################


sub getEntry {
    local ($id) = @_;

    $expEntry{$id}[$expNentries{$id}-1];
}


############################################################
# &addEntry ($lineType, $value)
#
# Adds a new entry to experiment file structure
# 
############################################################


sub addEntry {
    local ($id,$value) = @_;
    local ($n);

    if ($value !~ /\n$/) { $value .= "\n"; }

    $n = $expNentries{$id};
    $expEntry{$id}[$n] = $value;
    $expNentries{$id}++;
    push(@expOrder,$id); #keep a track of order
}

############################################################
# &replaceEntry ($id, $value)
#
# Removes any existing entries of $id
# Adds a new entry.
# 
############################################################

sub replaceEntry {
    local ($id,$value) = @_;

    &killEntry($id);
    &addEntry($id,$value);
}

############################################################
# &killEntry ($id)
#
# Removes all existing entries of $id
# 
############################################################

sub killEntry {
    my ($id) = @_;
    my($x);

    # Erase any expfile entries for this $id
    undef ($expEntry{$id});

    # Set counter to one (zero)
    $expNentries{$id}=0;

    # Spice out all previous $id entries
    for ($x=0;$x<=$#expOrder;$x++) {
	if($expOrder[$x] eq $id) {
	    splice(@expOrder,$x,1);
	    $x--;
	}
    }
}


    
	     

############################################################
# &writeExpFile(FILEHANDLE)
#
# Writes the experiment file to the FILEHANDLE
# The order in which lines are added is in the original order.
# YOU SHOULD NOT add extra entries except through &addEntry($type,$value);
#
###################################################




sub writeExpFile {
    local($file) = @_ ;
    local($i);
    local($id);
    local(%count);
    local ($seq, $av_count);

    # always have ID line first
    if (defined($id = &getEntry("ID"))) {
	print $file "ID   $id";
    }

    for($i = 0; $i <= $#expOrder; $i++) {
	$id = $expOrder[$i];
	if ($id eq "ID") {
	    # Skip
	} elsif ($id eq "SQ") {
	    ############################################################
	    # special treatment for sequences
	    ############################################################
	    print $file "$id\n";
	    $seq = $expEntry{$id}[$count{$id}];
	    $seq =~ s/\s//g; # Kill white space
	    $len = length($seq);
	    for ($j = 0 ; $j < $len ; ) {
		print $file "    ";
		for($k=0;$k<6 && $j<$len;$k++,$j+=10) {
		    print $file " ",substr($seq, $j, 10);
		}
		print $file "\n";
	    }
	    print $file "//\n";
	} elsif ($id eq "AV") {
	    ############################################################
	    # different special treatment for quality measures
	    ############################################################
            $av_count = 1;
	    my @avVals = split(/\s+\s*/, $expEntry{$id}[$count{$id}]);
	    #foreach ( split(/\n/,$expEntry{$id}[$count{$id}]) ) {
	    while (@avVals) {
		my $avLine = join(" ", splice(@avVals, 0, 20));
                if ($av_count) {
                   print $file "$id   $avLine\n";
                   $av_count = 0;
                } else {
		    print $file "$id        $avLine\n";
                }
	    }
	} else {
	    ############################################################
	    # treat everything else as a multiline
	    ############################################################
	    foreach ( split(/\n/,$expEntry{$id}[$count{$id}]) ) {
		print $file "$id   $_\n";
	    }
	}
	$count{$id}++;
    }
}




############################################################
# ($left, $right) = &determineActiveRegion()
#
# Determine the active (unclipped) region as specified in an
# experiment file, previously read in by &readExpFile()
#
# This algorithm is borrowed from xgap.
#
############################################################

sub determineActiveRegion {
    local($SQ);
    local($SQlen);
    local($CS);
    local($CSfrom);
    local($CSto);
    local($QL);
    local($QR);
    local($SL);
    local($SR);
    local($l);
    local($r);

    # Determine left and right end cuts
    if (! defined($SQ = &getEntry("SQ"))) {
	warn "determineActiveRegion - No sequence\n";
	$SQlen = 4096;
    } else {
	$SQlen = length($SQ);
    }
    if (! defined($SL = &getEntry("SL"))) { $SL = 0; }
    if (! defined($SR = &getEntry("SR"))) { $SR = $SQlen+1; }
    if (! defined($QL = &getEntry("QL"))) { $QL = 0; }
    if (! defined($QR = &getEntry("QR"))) { $QR = $SQlen+1; }

    #/* find where the good bits are */
    if ($SL>$QL) { $l = $SL; } else { $l = $QL; }
    if ($SR<$QR) { $r = $SR; } else { $r = $QR; }

    if ( defined($CS = &getEntry("CS"))) {

        ($CSfrom,$CSto) = $CS =~ /^(\d+)..(\d+)$/;

        #/*
        # * Six posibilities:
        # *           l                r
        # * 1  -----  >                <
        # * 2  -------|------->        <
        # * 3  -------|<---------------|-------
        # * 4         >    <--------   |
        # * 5         >    <-----------|-------
        # * 6         >                <  -----
        # *
        # * key: [-] CS [>] new left c.o [<] new right c.o [|] orig c.o
        # *
        # * Coded for understanding not efficiency!!!
        # */
        if ($CSfrom <= $l && $CSto <= $l) { #/*1*/
            $l = $l; #/* null statement */
        } elsif ($CSfrom <= $l+1 && $CSto < $r) { #/*2*/
            $l = $CSto;
        } elsif ($CSfrom <= $l+1 && $CSto >= $r) { #/*3*/
            $r = $l+1;
        } elsif ($CSfrom < $r && $CSto < $r) { #/*4*/
            $r = $CSfrom;
        } elsif ($CSfrom < $r && $CSto >= $r) { #/*5*/
            $r = $CSfrom;
        } elsif ($CSfrom >= $r && $CSto >= $r) { #/*6*/
            $l = $l; #/* null statement */
	}
    }
    if ($l>$r) { $l = $r-1; } #/* just a quick consistency check */

    ($l+1 , $r-1);
}

############################################################
############################################################
###  Routines to determine uncovered regions in sequences, contigs
############################################################
############################################################

############################################################
# &dumpcover(@uncovered)
#
# Dumps out a coverage array to STDERR
#
############################################################

sub dumpcover {
    local ($i);
    local (@uncovered) = @_;

    warn $#uncovered+1," uncovered regions\n";
    for ($i = 0; $i <= $#uncovered; $i++) {
	warn "  $uncovered[$i]->[0] $uncovered[$i]->[1]\n";
    }
}


############################################################
# &initcover($from, $to)
#
# Initialised array @uncovered with $from and $to. This is
# normally called with $from=1 and $to=length($seq).
#
############################################################
sub initcover {
    local ( $c1 , $c2 ) = @_;

    @uncovered = ();

    if ($c1 > $c2) { ($c1, $c2) = ($c2, $c1); }

    $uncovered[0] = [$c1, $c2];

}

############################################################
# &cover($from, $to)
#
# Mask out region $from..$to from @uncovered array.
#
############################################################
sub cover {
    local ( $c1 , $c2 ) = @_;
    local($i);

    if ($c1 > $c2) { ($c1, $c2) = ($c2, $c1); }

    for ($i = 0; $i <= $#uncovered; $i++) {
	# Does $c1..$c2 cover $uncovered[$i]->[0]..$uncovered[$i]->[1] ???
	# Consider all 6 actions
	if ($c2 < $uncovered[$i]->[0]) {
	    last;
	} elsif ($c1 <= $uncovered[$i]->[0] && $c2 < $uncovered[$i]->[1]) {
	    $uncovered[$i]->[0] = $c2+1;
	} elsif ($c1 <= $uncovered[$i]->[0]) {
	    splice(@uncovered,$i,1); $i--;
	} elsif ($c2 < $uncovered[$i]->[1]) {
	    splice(@uncovered,$i,0,[$uncovered[$i]->[0],$uncovered[$i]->[1]]);
	    $uncovered[$i]->[1] = $c1-1;
	    $i++;
	    $uncovered[$i]->[0] = $c2+1;
	} elsif ($c1 <= $uncovered[$i]->[1]) {
	    $uncovered[$i]->[1] = $c1-1;
	} else {
	    # no action
	}
    }

}

############################################################
# ($left, $right) = &myDetermineActiveRegion($mode)
#
# This is a rewrite of &determineActiveRegion() using &cover()
#
# Mode:
#     -c  - mask cloning vector (CS)
#     -C  - mask everything but the cloning vector
#     -ql - mask poor quality regions (1..QL)
#     -qr - mask poor quality regions (QR..END)
#     -sl - mask sequencing vector (1..SL)
#     -sr - mask sequencing vector (SR..END)
#     -xl - mask using quality revert Squirrels final clip points
#     -xr - mask using quality revert Squirrels final clip points
#
############################################################

sub myDetermineActiveRegion {
    local($SQ);
    local($SQlen);
    local($CS);
    local($CSfrom);
    local($CSto);
    local($QL);
    local($QR);
    local($exQL);
    local($exQR);
    local($SL);
    local($SR);
    local(@overall);
    local(@uncovered);
    local($switches) = @_;

    # Determine left and right end cuts
    if (! defined($SQ = &getEntry("SQ"))) {
	warn "determineActiveRegion - No sequence\n";
	$SQlen = 4096;
    } else {
	$SQlen = length($SQ);
    }


    # determine active region
    &initcover(1,$SQlen);
    if ($switches =~ /-sl/) { 
	if (defined($SL = &getEntry("SL"))) { &cover(0,$SL); }
    }
    if ($switches =~ /-sr/) { 
	if (defined($SR = &getEntry("SR"))) { &cover($SR,$SQlen+1); }
    }
    if ($switches =~ /-ql/) { 
	if (defined($QL = &getEntry("QL"))) { &cover(0,$QL); }
    }
    if ($switches =~ /-qr/) {
	if (defined($QR = &getEntry("QR"))) { &cover($QR,$SQlen+1); }
    }
    if ($switches =~ /-c/) {
	if ( defined($CS = &getEntry("CS"))) {
	    ($CSfrom,$CSto) = $CS =~ /^(\d+)..(\d+)$/;
	    &cover($CSfrom,$CSto);
	}
    }
    if ($switches =~ /-C/) {
	if ( defined($CS = &getEntry("CS"))) {
	    ($CSfrom,$CSto) = $CS =~ /^(\d+)..(\d+)$/;
	    &cover(0,$CSfrom-1);
	    &cover($CSto+1,$SQlen+1);
	}
    }

    ############################################################
    # Precheck
    if ($switches =~ /-x/) {
	# We should ignore overall clips!
	$overall = 0;
	foreach ( @{$expEntry{"CC"}} ) {
	    if (/Set final clip points for gap/) {
		$overall = 1;
		last;
	    }
	}
    }
    ############################################################

    if ($switches =~ /-xl/) {
	if (defined($exQL = $expEntry{"QL"}[$expNentries{"QL"}-1-$overall])) { &cover(0,$exQL); }
    }
    if ($switches =~ /-xr/) {
	if (defined($exQR = $expEntry{"QR"}[$expNentries{"QR"}-1-$overall])) { &cover($exQR,$SQlen+1); }
    }


    # Return results
    if ($#uncovered >= 0) {
	return ($uncovered[0]->[0],$uncovered[0]->[1]);
    } else {
	return (1,0);
    }
}

1;

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

