#! /usr/local/bin/perl
#
# Usage: parse_notes.pl <full assembly.notes filename+path> <near_complete>
#
# Reads through an assembly.notes file and makes entries in Oracle of 
# the following:
# project_phrap entries: one per detected phrap run, incl insert size and q20
# phrap_markers: one entry for each of the latest set of sts marker hits/misses
# clone_overlaps: updates and inserts overlap entries from cross_match hits 
# (see notes in the code for the rules)
# If the optional second commandline argument is provided, all 
# overlaps not yet examined (confirmed or refuted) are marked as refuted 
# (This switch should only be included when the project is close to completion, 
# and any overlap not detected here is probably wrong).
#
# SJM, 28.4.00
#
# Configuration ----------------------------------------------------------

# $Id: parse_notes.pl,v 1.19 2002/05/30 13:48:50 rmd Exp $

BEGIN {
	$ENV{'TWO_TASK'} = 'sids';
	$ENV{'ORACLE_HOME'} = '/usr/local/oracle';	# Oracle home
	push(@INC,'/usr/local/badger/bin');		# Path for WrapDBI.pm
}

use strict;
use WrapDBI;		# Sanger wrapper for the DBperl interface
# use Data::Dumper;	# For debugging only

my $dbi_user = 'assembly';	# DBI alias
my $debug = 0; # 1 to activate verbose reporting
my $q20_switch;
my $sts_switch;
my $overlap_switch;
my $checksanger_file; # Name of file to parse for HUMAN checksanger output
my $dbh;
my %phraprun; # Data on a single run
my @phrap_runs; # Array of phrap run hashes
my @marker_table; # Array of hashes for the latest marker results
my @ol_table; # Array of hashes for the latest cross_match results

# Location of logfile for error and warning messages
my $logfile = '/nfs/disk65/sm1/parse_notes.log'; 
# Don't report overlaps with these IDs
my $exclusions_file = '/usr/local/badger/bin/exclude_ids.txt'; 
my $ol_size = 2500; # Minimum size in bp before an overlap is recorded
my $pct_homol = 97; # Minimum % homology before an overlap is recorded
my $cs_threshold = 1400; # Checksanger matches recorded when more than this

# Main -------------------------------------------------------------------

my ($notesraw,$near_complete) = @ARGV;

my ($notespath,$notesfile) = $notesraw =~ /(.+\/)([^\/]+)$/;
my ($project) = $notesfile =~ /^(\S+).assembly.notes/;
chomp(my $timestamp = `date`);

unless ($project) {
	&write_log("ERROR: Could not find the projectname in path: ($notesraw)");
	die "ERROR: parse_notes could not find the projectname
path: ($notesraw)
This will not affect phrap2gap, but some data will not be stored in Oracle\n";
}

if ($near_complete) {
	&write_log("$timestamp Parsing $project (Near completion)");
} else {
	&write_log("$timestamp Parsing $project (Standard run)");
}

# get list of overlaps to exclude from Oracle
my @exclude_ids = &parse_excludes($exclusions_file);

open (NOTES,$notesraw) || die "Could not find $notesraw\n";
while (<NOTES>) {
	chop;
	unless ($_) {
		# Blank line - reset parsing switches
		$q20_switch = 0;
		$sts_switch = 0;
		$overlap_switch = 0;
	}
	if (/^Project: (\S+)$/) {
		unless ($1 eq $project) {
			print STDERR "ERROR: Project $1 not the same as $project detected in $notesfile\n";
			&write_log("ERROR: Project $1 not the same as $project detected in $notesfile");
		}
	}

	if (/^Started at \S+\s+(\S+)\s+(\d+) (\d+):(\d+):(\d+) \S+ (\d+)/) {
		# Start of a phrap run
		# Examples
		# 	Started at Tue Feb  1 19:53:54 GMT 2000
		# 	Started at Tue Feb 16 11:43:42 GMT 1999
		# Reformat date string to match the Oracle format of 'Mon-DD-YYYY:HH24:MI:SS'
		print STDERR "Date line: $_\n" if $debug;
		$phraprun{'Run_date'} = "$1-$2-$6:$3:$4:$5";
	}
	
	if (/^Phred 20 coverage estimate for /) {
		# detects start of Phred 20 (Q20) results
		$q20_switch = 1;
		$sts_switch = 0;
		$overlap_switch = 0;
	} elsif ($_ eq '------------ --------------- ----------- ----------') {
		# detects start of table of sts results (one per phrap2gap run with the -sts switch)
		# only record the latest set of results
		undef @marker_table;
		$q20_switch = 0;
		$sts_switch = 1;
		$overlap_switch = 0;
	} elsif ($_ eq 'Clone           Length Matches Length Identity') {
		# detects the start of the overlaps table
		# only record the latest set of results
		undef @ol_table;
		$q20_switch = 0;
		$sts_switch = 0;
		$overlap_switch = 1;
	} elsif (/^cafcheck2: Please see file (\S+) for the full match list/) {
		# gets the name of the checksanger file
		# (matches more than once but only retains the final filename)
		$checksanger_file = $1;
		print STDERR "Checksanger file is now ($1)\n" if $debug;
	}

	if ($q20_switch) {
		print STDERR "Q20 line: $_\n" if $debug;
		&parse_phraprun($_,$project);
	} elsif ($sts_switch) {
		print STDERR "Marker line: $_\n" if $debug;
		&parse_markerline($_);
	} elsif ($overlap_switch) {
		print STDERR "Overlap line: $_\n" if $debug;
		&parse_overlapline($_, $ol_size, $pct_homol);
	}
}
close(NOTES);

# get checksanger results
my @csresults = &parse_checksanger($notespath,$checksanger_file);

# For debug only
# print Data::Dumper::Dumper(\@marker_table), "\n";
# print Data::Dumper::Dumper(\@ol_table), "\n";
# print Data::Dumper::Dumper(\@phrap_runs), "\n";

eval {
    &write_log("Writing $project");
    $dbh = WrapDBI->connect($dbi_user, {RaiseError => 1});
    my $get_clonename_sth = $dbh->prepare(qq[select clonename from clone_project where projectname = ?]);
    my $phrap_ins_sth = $dbh->prepare(qq[insert into project_phrap (projectname,iscurrent,phrap_date,est_ins_size_kb,q20_val_1kb)values (?,0,to_date(?, 'Mon-DD-YYYY:HH24:MI:SS'),?,?)]);
    my $iscurrent_sth1 = $dbh->prepare(qq[update project_phrap set iscurrent = 0 where projectname = ?]);
    my $iscurrent_sth2 = $dbh->prepare(qq[update project_phrap set iscurrent = 1 where projectname = ? and phrap_date = (select max(phrap_date) from project_phrap where projectname = ?)]);
    my $mar_ext_sth = $dbh->prepare(qq[select marker_name from phrap_markers where projectname = ?]);
    my $mar_del_sth = $dbh->prepare(qq[delete from phrap_markers where projectname = ?]);
    my $mar_ins_sth = $dbh->prepare(qq[insert into phrap_markers (projectname,marker_name,marker_pos,on_contig,expected) values (?,?,?,?,?)]);
    my $ol_ext_sth = $dbh->prepare(qq[select clonename,overlaps_clone,confirmation,remark from clone_overlaps where clonename = ? or overlaps_clone = ?]);
    my $ol_upd_sth = $dbh->prepare(qq[update clone_overlaps set size_bp = ? where clonename = ? and overlaps_clone = ?]);
    my $ol_upd_sth2 = $dbh->prepare(qq[update clone_overlaps set remark = ? where clonename = ? and overlaps_clone = ?]);
    my $ol_ins_sth = $dbh->prepare(qq[insert into clone_overlaps (clonename,overlaps_clone,confirmation,size_bp,remark)  values(?,?,?,?,?)]);
    my $ol_unc_sth = $dbh->prepare(qq[update clone_overlaps set confirmation = 0, remark = ? where confirmation is null and (clonename = ? or overlaps_clone = ?)]);
    
    # find the clone name for this project
    $get_clonename_sth->execute($project);
    my $clonename = ($get_clonename_sth->fetchrow_array)[0];
    
    if ($phrap_runs[0]{'Project'}) {
	my $made_inserts = 0;
	foreach my $pref (@phrap_runs) {
	    # Try to insert any new phrap runs - ignore duplicate insert errors
	    eval {
		if ($pref->{'Project'} && $pref->{'Run_date'} && $pref->{'Insert_size'} && $pref->{'Q20_val'}) {
		    print STDERR "Inserting Phrap run for $pref->{'Project'} on $pref->{'Run_date'}\n" if $debug;
		    $phrap_ins_sth->execute($pref->{'Project'},$pref->{'Run_date'},$pref->{'Insert_size'},$pref->{'Q20_val'});
		    $made_inserts = 1;
		} else {
		    print STDERR "Incomplete Phrap run record for $pref->{'Project'}\n" if $debug;
		    &write_log("Incomplete Phrap run record for $pref->{'Project'}");
		}
	    };
	    if ($@) {
		if ($DBI::err ==1) {
		    print STDERR "Tried to insert a duplicate phrap run record - ignored\n" if $debug;
		} elsif ($DBI::err) {
		    &write_log($DBI::errstr);
		    die "ERROR: $DBI::errstr";
		} else {
		    &write_log($@);
		    die "ERROR: $@\n";
		}
	    }
	}
	if ($made_inserts) {
	    # Only do these updates if at least one new insert has been made
	    # Mark all project_phrap records as non-current
	    $iscurrent_sth1->execute($project);
	    # Mark latest project_phrap record as current
	    $iscurrent_sth2->execute($project, $project);
	}
    }
    
    unless($clonename) {
	# Marker and overlap information is associated with the clone, not the project.
	# No clonename means this project is not linked to a clone in Oracle, hence these types of
	# data cannot be inserted into Oracle and this program exits here.
	my $err_str = "Oracle has no clone record associated with project $project.
Overlap and marker information cannot therefore be inserted into Oracle
Please enter the associated clone record into Oracle (using Gull) or contact bug.\n";
	die $err_str;
    }
    
    if ($marker_table[0]{'Project'}) {
	# Markers were found in the assembly.notes file
	# Duplicates will cause parsing to stop at that point in the table - this is okay
	# (duplicates are a phrap2gap bug and occur at the end of the table)
	
	eval {
	    $mar_ext_sth->execute($project);
	    if (($mar_ext_sth->fetchrow_array)[0]) {
		# Markers exist for this project, and a set have been found
		# in the assembly.notes file - delete and replace them
		$mar_del_sth->execute($project);
		print STDERR "Deleting markers for $project\n" if $debug;
	    }
	    foreach my $mref (@marker_table) {
		print STDERR "Inserting Phrap marker $mref->{'Marker'} for $mref->{'Project'}\n" if $debug;
		$mar_ins_sth->execute($mref->{'Project'},$mref->{'Marker'},$mref->{'Matchpos'},$mref->{'Contig'},$mref->{'Expected'});
	    }
	};
	if ($@) {
	    if ($DBI::err ==1) {
		print STDERR "Tried to insert a duplicate marker - ignored\n" if $debug;
	    } elsif ($DBI::err) {
		&write_log($DBI::errstr);
		die "ERROR: $DBI::errstr";
	    } else {
		&write_log($@);
		die "ERROR: $@\n";
	    }
	}
    }
    
    if ($ol_table[0]{'Project'}) {
	# A table of cross_match results was found in the assembly.notes file
	# What happens now:
	# Each overlap is case-corrected to compensate for checksanger's 
	#    uppercasing (this is for backward compatibility - checksanger 
	#    should not uppercase now)
	# Overlaps are then rejected if a match is detected to itself 
	#    or to its EMBL entry (if this sequence is matched against 
	#    itself but the names are not similar, this self-match will 
	#    not be detected)
	# Overlaps are rejected if they are with a set of identifiers 
	#    in a file (vector sequences, etc)
	# If an overlap is found which is already known to Oracle:
	# - if it has been confirmed or refuted, leave it alone
	# - if it has not been examined, update the overlap size
	# If an overlap is found which has not previously been found:
	# - an overlap entry is created, with a remark ('Found by cross_match') 
	#   but with a null confirmation (unexamined by finisher)
	# - Overlaps not found by this program will not be touched in the 
	#   default mode BUT if run with the near_complete flag set, all 
	#   unconfirmed overlaps will be refuted - THIS MODE MAY BE DISABLED
	# - Overlaps are only represented in one direction, i.e, if 
	#   A overlaps B is in the database, B overlaps A will not be inserted, 
	#   but the already extant overlap will be confirmed

	foreach my $olref (@ol_table) {
	    # Change case of matched name if necessary (retained for backwards compatibility)
	    $olref->{'Matches_clone'} = &format_checksanger($olref->{'Matches_clone'});
	}
	
	foreach my $olref (@ol_table) {
	    # Reject any overlap between this clone and an ID on the exclusions list
	    foreach my $exclude (@exclude_ids) {
		if ($olref->{'Matches_clone'} eq $exclude) {
		    print STDERR "Rejected overlap with $olref->{'Matches_clone'}: ID is on the exclusion list\n" if $debug;
		    $olref->{'reject'} = 1;
		    next;
		}
	    }
	}
	
	foreach my $olref (@ol_table) {
	    # Flag duplicate overlaps (ie, matches found with the EMBL entry and the same Sanger project)
	    my ($curr_fullname, $curr_prefix, $curr_name) = $olref->{'Matches_clone'} =~ /(([A-Za-z]+)(\S+))/;
	    my $curr_matchlen = $olref->{'Match_len'};
	    foreach my $oref (@ol_table) {
		my ($fullname, $prefix, $name) = $oref->{'Matches_clone'} =~ /(([A-Za-z]+)(\S+))/;
		if  ($curr_name eq $name && $curr_fullname ne $fullname) {
		    # not me - something similar
		    if ($curr_matchlen == $oref->{'Match_len'} && $curr_prefix =~ /HS|HSDJ/i) {
			# EMBL prefix but same overlap size - reject
			$olref->{'reject'} = 1;
			print STDERR "Rejected overlap with $curr_fullname: duplicate of $fullname\n" if $debug;
			next;
		    }
		} elsif ($curr_fullname eq $clonename && ! $olref->{'reject'}) {
		    $olref->{'reject'} = 1;
		    print STDERR "Rejected overlap with $curr_fullname for being the clone ($clonename)\n" if $debug;
		    next;
		}
	    }
	}
	
	# find existing overlaps (in the clone_overlaps table)
	$ol_ext_sth->execute($clonename,$clonename);
	my $ol_arrayref = $ol_ext_sth->fetchall_arrayref;
	
	foreach my $olref (@ol_table) {
	    next if $olref->{'reject'};
	    my $found_overlap = 0;
	    foreach(@{$ol_arrayref}) {
		# Checks new overlap against all current entries, updating if necessary.
		# values returned in this order:
		#   clonename,overlaps_clone,confirmation,remark
		if ($_->[0] eq $clonename && $_->[1] eq $olref->{'Matches_clone'}) {

		    # overlap found - in forward direction
		    $found_overlap = 1;

		    if ($_->[2] eq '') {
			# update size_bp and remark if confirmation is null (unexamined overlap)

			print STDERR "Updating overlap between $clonename and $olref->{'Matches_clone'}\n" if $debug;

			$ol_upd_sth->execute($olref->{'Match_len'},$clonename,$olref->{'Matches_clone'});

			unless ($_->[3] =~ /Found by cross_match/i) {
			    # Update comment unless it has been done already
			    $ol_upd_sth2->execute($_->[3].'(also found by cross_match)',$clonename,$olref->{'Matches_clone'});
			}
		    }
		} elsif ($_->[1] eq $clonename && $_->[0] eq $olref->{'Matches_clone'}) {

		    # overlap found - in reverse direction

		    $found_overlap = 1;

		    if ($_->[2] eq '') {
			# update size_bp and remark if confirmation is null (unexamined overlap)

			print STDERR "Updating reverse overlap between $clonename and $olref->{'Matches_clone'}\n" if $debug;

			$ol_upd_sth->execute($olref->{'Match_len'},$olref->{'Matches_clone'},$clonename);

			unless ($_->[3] =~ /Found by cross_match/i) {
			    # Update comment unless it has been done already
			    $ol_upd_sth2->execute($_->[3].'(also found by cross_match)',$olref->{'Matches_clone'},$clonename);
			}
		    }
		}
	    }

#	    unless ($found_overlap) {
#		# overlap does not already exist
#		foreach (@csresults) {
#		    if ($_ eq $olref->{'Matches_clone'}) {
#			# overlap exceeds threshold for contiguous match
#			print STDERR "Inserting overlap between $clonename and $olref->{'Matches_clone'}\n" if $debug;
#			$ol_ins_sth->execute($clonename,$olref->{'Matches_clone'},'',$olref->{'Match_len'},'Found by cross_match');
#			last;
#		    }
#		}
#	    }

	}
    }
    
    # mark unexamined overlaps as unconfirmed if this script is running on a
    # project close to completion, where all overlaps should be known
    $ol_unc_sth->execute('Automatically refuted: not found by cross_match',$clonename,$clonename) if $near_complete;
};

if ($@) {
    if ($DBI::err) {
		print "ERROR: $DBI::errstr\n";
		&write_log($DBI::errstr);
    }
	print "ERROR: $@\n";
	&write_log($@);
	if (defined $dbh) {$dbh->disconnect};
} else {
	$dbh->disconnect;
}

&write_log("Project $project complete\n");

# Subroutines ------------------------------------------------------------

sub parse_phraprun {
# Fills the @phrap_runs data structure with the details of each run
# Extracts Q20 and insert size from Phred paragraph
	my ($q20_line,$projectname) = @_;
	$phraprun{'Project'} = $projectname;
	if ($q20_line =~ /^Estimated insert size = (\S+) kb$/) {
		$phraprun{'Insert_size'} = $1;
	} elsif ($q20_line =~ /^1kb contigs (\S+)$/) {
		# Last line of the Q20 - write the phrap record into the data structure
		$phraprun{'Q20_val'} = $1;
		push @phrap_runs,{%phraprun};
		undef %phraprun;
	}
}

sub parse_overlapline {
# Extracts overlaps from the cross_match table
# If the line does not conform to the table structure, assume table end and stop parsing
#Clone           Length Matches Length Identity
	my ($ol_line,$ol_size,$pct_homol) = @_;
	my %olline;
	$ol_line =~ s/ +$//;
	my ($clone,$tot_len,$num_matches,$mean_len,$mean_id) = split(/\s+/,$ol_line,5);
	next if $clone eq 'Clone';
	if ($mean_id) {
		if ($mean_id =~/ /) {
			# line too long - assume we arent in the table anymore
			$overlap_switch = 0;
			next;
		}
	} else {
		# line too short - assume we arent in the table anymore
		$overlap_switch = 0;
		next;
	}
	$mean_id =~ s/%//;
	if (($tot_len >= $ol_size) && ($mean_id >= $pct_homol)) {
		# Overlap entry equals or exceeds the size and % identity thresholds
		$olline{'Project'} = $project;
		$olline{'Matches_clone'} = $clone;
		$olline{'Match_len'} = $tot_len;
		push @ol_table,{%olline};
	}
}

sub parse_markerline {
# Extracts marker matches and misses from the STS table
# If the line does not conform to the table structure, assume table end and stop parsing
#marker       match posistion  in contig   expected 
	my ($marker_line) = @_;
	$marker_line =~ s/no match/no_match/g;
	$marker_line =~ s/ +$//;
	my %markerline;
	my ($marker,$matchpos,$contig,$expected) = split(/\s+/,$marker_line,4);
	next if $marker eq '------------';
	if ($expected) {
		if ($expected =~/ /) {
			# line too long - assume we arent in the table anymore
			$sts_switch = 0;
			next;
		}
	} else {
		# line too short - assume we arent in the table anymore
		$sts_switch = 0;
		next;
	}

	if ($expected eq 'YES') {
		$expected = 1;
	} else {
		$expected = 0;
	}	
	if ($contig eq '******') {$contig = ''};

	$markerline{'Project'} = $project;
	$markerline{'Marker'} = $marker;
	$markerline{'Matchpos'} = $matchpos;
	$markerline{'Contig'} = $contig;
	$markerline{'Expected'} = $expected;
	push @marker_table,{%markerline};
}

sub write_log {
# Writes errors and warnings to a logfile
# Prints to STDERR if it can't write to the logfile
	my ($error) = @_;
#	if (open LOG, ">>$logfile") {
#		print LOG "$error\n";
#		close LOG;
#	} else {
#		print STDERR "ERROR: Can't write to logfile $logfile\n";
		print "$error\n";
#	}
}

sub parse_checksanger {
# Parses the checksanger file for matches greater than the checksanger threshold value
	my ($projpath,$cs_name) = @_;
	my @csresults;
	my $csfilename = "$projpath$cs_name";
	if (open CSFILE,$csfilename) {
		while (<CSFILE>) {
			my @csline = split;
			if ($csline[3] > $cs_threshold) {
			my $clonename = &format_checksanger($csline[0]);
				push(@csresults,$clonename);
				print STDERR "Found significant overlap with $clonename (unformatted: $csline[0]) - $csline[3] exceeds $cs_threshold\n" if $debug;
			}
		}
		close(CSFILE);
	} else {
		print STDERR "Could not open $csfilename\n" if $debug;
	}
	return @csresults;
}

sub parse_excludes {
# Parses the exclude IDs file and returns an array of the IDs to exclude
	my ($exclude_file) = @_;
	my @excludes;
	if (open EXFILE,$exclude_file) {
		while (<EXFILE>) {
			chop;
			next unless $_;
			next if /$\s*#/;
			my ($ex_id,$rest) = split(/\s+/,$_,2);
			print STDERR "Exclude overlaps with ID: $ex_id\n" if $debug;
			push(@excludes,$ex_id);
		}
		close(EXFILE);
	} else {
		print STDERR "Could not open $exclude_file\n" if $debug;
	}
	return @excludes;
}

sub format_checksanger {
# Stupidly, checksanger uppercases all clone names!
# This attempts to rectify the situation
# (checksanger now fixed - does not uppercase names, but this code still needed to parse old output)
	my ($cs_clone) = @_;
	$cs_clone =~ s/^ST(\d)/St$1/o;
	$cs_clone =~ s/^B(\d)/b$1/o;
	$cs_clone =~ s/^BA(\d)/bA$1/o;
	$cs_clone =~ s/^BB(\d)/bB$1/o;
	$cs_clone =~ s/^BG(\d)/bG$1/o;
	$cs_clone =~ s/^BK(\d)/bK$1/o;
	$cs_clone =~ s/^BM(\d)/bM$1/o;
	$cs_clone =~ s/^C(\d)/c$1/o;
	$cs_clone =~ s/^CA(\d)/cA$1/o;
	$cs_clone =~ s/^CB(\d)/cB$1/o;
	$cs_clone =~ s/^CE(\d)/cE$1/o;
	$cs_clone =~ s/^CI(\d)/cI$1/o;
	$cs_clone =~ s/^CL(\d)/cL$1/o;
	$cs_clone =~ s/^CN(\d)/cN$1/o;
	$cs_clone =~ s/^CU(\d)/cU$1/o;
	$cs_clone =~ s/^CV(\d)/cV$1/o;
	$cs_clone =~ s/^DA(\d)/dA$1/o;
	$cs_clone =~ s/^DJ(\d)/dJ$1/o;
	$cs_clone =~ s/^DM(\d)/dM$1/o;
	$cs_clone =~ s/^FA(\d)/fA$1/o;
	$cs_clone =~ s/^FF(\d)/fF$1/o;
	$cs_clone =~ s/^P(\d)/p$1/o;
	return $cs_clone;
}

__END__

=head1 Log
 # $Log: parse_notes.pl,v $
 # Revision 1.19  2002/05/30  13:48:50  rmd
 # Turned dubugging off!
 #
 # Revision 1.18  2002/05/30  13:47:21  rmd
 # Stopped putting in overlaps found only by cross_match as there were far too
 # many of them.
 # Added log.
 #
