#!/usr/local/bin/perl5.004
#######################################################################
# 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: proj_depth.pl,v 1.6 2007/01/17 16:20:14 ja1 Exp $	

use strict;
use WrapDBI;
use Data::Dumper;
use SDTable;

my $debugging = shift || 0;

warn "$debugging\n" if $debugging;

sub numerically { $a <=> $b; }

Search(100.);

sub Search {
    
    my ($Q20_thresh) = @_;
    unless ($Q20_thresh) { return; }
    unless ($Q20_thresh =~ /\d+/) { return; }


    my $db = "reports";
    my $dbh = WrapDBI->connect("$db", {RaiseError => 1});

# statuses are 'Shotgun','Assembly start','Half shotgun complete', 'Shotgun complete', 'Selected for top-up'    
#    my $statement = "
#select distinct a.clonename, c.projectname, g.teamname, i.description, d.description
#from clone a, clone_project b, project_status c, projectstatusdict d, clone_priority e, cloneprioritydict f, project_owner g, projectownerdict i
#where a.clonename = b.clonename
#and b.projectname = c.projectname
#and c.status = d.id_dict
#and a.clonename = e.clonename
#and e.priority = f.id_dict
#and b.projectname = g.projectname
#and g.owner_type = i.id_dict
#and a.sequenced_by = 5
#and e.is_current = 1
#and g.iscurrent = 1
#and c.iscurrent = 1
#and c.status in (13,16,30,15,33)
#order by d.description, g.teamname, i.description, c.projectname
#";

my $statement = "
select distinct a.clonename, c.projectname, g.teamname, i.description,
d.description, a.sequenced_by
from clone a
, clone_project b
, project_status c
, projectstatusdict d
, clone_priority e
, cloneprioritydict f
, project_owner g
, projectownerdict i
where a.clonename = b.clonename
and b.projectname = c.projectname
and c.status = d.id_dict
and a.clonename = e.clonename
and e.priority = f.id_dict
and b.projectname = g.projectname
and g.owner_type = i.id_dict
and e.is_current = 1
and g.iscurrent = 1
and c.iscurrent = 1
and c.status in (13,16,30,15,33)
order by d.description, g.teamname, i.description, c.projectname
";
    
#    warn "$statement\n";
    warn scalar localtime(), "\n";
    
    my $sth1 = $dbh->prepare($statement) || die "ERROR:$DBI::errstr";
    $sth1->execute || die "ERROR:$DBI::errstr";
    
    my ($clonename, $projectname, $teamname, $description, $status, $by);
    my $hash = {};
    my $count = 0;
    while ( ($clonename, $projectname, $teamname, $description, $status, $by) = $sth1->fetchrow() ){
	next unless $by == 5;
	$count++;
#	warn $count;
	my @array = ($clonename, $projectname, $teamname, $description, $status);
	$hash->{$count} = \@array;
	last if $debugging and $count >= 50;
    }
    $sth1->finish();

    warn scalar localtime(), ": $count clones found\n";


    
    ############################################
    # Get Q20 figures from project_phrap tables
    ############################################
    my $statement4 = "
select a.EST_INS_SIZE_KB, a.Q20_VAL_1KB
from project_phrap a
where a.PROJECTNAME = ?
and a.iscurrent = 1
";
    my ($EST_INS_SIZE_KB, $Q20_VAL_1KB, $q_proj);
    my $sth4 = $dbh->prepare($statement4) || die "ERROR:$DBI::errstr";
    foreach my $n ( sort numerically ( keys(%$hash) ) ){

	$q_proj = $hash->{$n}->[1];
	$sth4->bind_param(1, $q_proj);
	$sth4->execute || die "ERROR:$DBI::errstr";
	while ( ($EST_INS_SIZE_KB, $Q20_VAL_1KB) = $sth4->fetchrow() ){
 	    $hash->{$n}->[5] = $EST_INS_SIZE_KB;
	    $hash->{$n}->[6] = $Q20_VAL_1KB;
	}
    }
    $sth4->finish();
    warn scalar localtime(), ": Got Q20 figures\n";
    #######################################
    # Get some project_summary information
    #######################################

    my $statement5 = "
select num_lge_contigs, lge_ctg_len, num_reads
from project_summary 
where projectname = ? and summarydate in
(select max(summarydate) from project_summary a 
where a.projectname = projectname)
";
#select a.NUM_LGE_CONTIGS, LGE_CTG_LEN, NUM_READS
#from project_summary a
#where a.PROJECTNAME = ?
#and (a.projectname, a.summarydate) in (
#    select b.projectname, max(b.summarydate)
#    from project_summary b
#    group by b.projectname)
#";
    #######################################
    # Get restriction digest size info
    #######################################
    my $statement6 = "
select min(a.insert_size_bp), round(avg(a.insert_size_bp)), max(a.insert_size_bp), count(*)
from rdgel_lane_image a, 
     rdgel_lane b, 
     rdrequest c 
where a.id_rdgel = b.id_rdgel 
  and a.lane = b.lane 
  and b.id_rdrequest = c.id_rdrequest
  and a.isgood = 1
  and c.clonename = ?
";
    my $statement7 = "
select b.speciesname || ' ' || b.chromosome 
from clone a, chromosomedict b
where b.id_dict = a.chromosome
and a.clonename = ?
";
    my ($NUM_LGE_CONTIGS, $LGE_CTG_LEN, $NUM_READS, $s_proj, $s_clone);
    my $sth5 = $dbh->prepare($statement5) || die "ERROR:$DBI::errstr";
    my $sth6 = $dbh->prepare($statement6) || die "ERROR:$DBI::errstr";
    my $sth7 = $dbh->prepare($statement7) || die "ERROR:$DBI::errstr";
    foreach my $n ( sort numerically ( keys(%$hash) ) ){

        $s_proj = $hash->{$n}->[1];
        $sth5->bind_param(1, $s_proj);
        $sth5->execute || die "ERROR:$DBI::errstr";
        while ( ($NUM_LGE_CONTIGS, $LGE_CTG_LEN, $NUM_READS) = $sth5->fetchrow()){
            $hash->{$n}->[7] = $NUM_LGE_CONTIGS;
            $hash->{$n}->[8] = $LGE_CTG_LEN;
            $hash->{$n}->[9] = $NUM_READS;
        }
        $s_clone = $hash->{$n}->[0];
        $sth6->bind_param(1, $s_clone);
        $sth6->execute || die "ERROR:$DBI::errstr";
	my ($minsize, $meansize, $maxsize, $numdigests) = $sth6->fetchrow();
	$hash->{$n}->[10] = sprintf("%d-%d-%d [%d]",
				    defined $minsize ? int($minsize/1000 + 0.5) : 0,
				    defined $meansize? int($meansize/1000 + 0.5) : 0,
				    defined $maxsize ? int($maxsize/1000 + 0.5) : 0,
				    defined $numdigests ? $numdigests : 0
				    );
        $sth7->bind_param(1, $s_clone);
        $sth7->execute || die "ERROR:$DBI::errstr";
	my ($chrom) =  $sth7->fetchrow();
	$hash->{$n}->[11] = $chrom;
    }
    $sth5->finish();
    $sth6->finish();
    $sth7->finish();

    warn scalar localtime(), ": Got project summary and digest info\n";


    my $table2 = SDTable->new('format' => 'gaudy');
    $table2->set('titles' => ['Team','Q20 Depth','Status','Project','Chrom','Size','Reads','Contigs',
			      'Digest (min-avg-max [n])']);
    
    my $new_count = 0;
    my $statement8 = "
select od.ONLINE_PATH
from ONLINE_DATA od, project p
where p.id_online = od.id_online
and od.IS_AVAILABLE = 1
and p.PROJECTNAME = ?
";
    my $sth8 = $dbh->prepare($statement8) || die "ERROR:$DBI::errstr";
    foreach my $n ( sort numerically ( keys(%$hash) ) ){
	if (defined $hash->{$n}->[6] and $hash->{$n}->[6] < $Q20_thresh){
	    $hash->{$n}->[7] .= " (1k)";
	    $hash->{$n}->[9] .= " (1k)";
	    my $s_proj = $hash->{$n}->[1];
	    $sth8->bind_param(1, $s_proj);
	    $sth8->execute || die "ERROR:$DBI::errstr";
	    my ($projpath) = $sth8->fetchrow();
	    if (defined $projpath) {
		my $notes = "$projpath/$s_proj.assembly.notes";
		if (open F, $notes) {
		    my @lines = <F>;
		    close F;
		    warn "$notes\n" if $debugging;
		    my @temp = grep /^(\|\s+Number of (contig|reading)s|Estimated insert size)/, @lines;
		    my ($contigs1k, $reads1k) = ($hash->{$n}->[7], $hash->{$n}->[9]);
		    my ($contigs2k, $reads2k) = (0,0);
		    warn @temp if $debugging;
		    my $i = $#temp;
		    my $found = 0;
		    my $pattern = sprintf "Estimated insert size = %.1f kb", $hash->{$n}->[5];
		    while (not $found and $i >= 0) {
			while ($i >= 0 and $temp[$i] !~ /$pattern/) {
			    $i--;
			}
			if ($i >= 0) {
			    foreach my $j ($i+1..$#temp) {
				last if $temp[$j] =~ /Estimated insert size = /;
				next unless $temp[$j] =~ /Number of contigs >2000\s+=\s+(\d+)/;
				$contigs2k = $1;
				if ($temp[$j+1] =~ /Number of readings\s+=\s+(\d+)/) {
				    $reads2k = $1;
				} 
				last;
			    }
			}
			if ($contigs2k) {
			    $hash->{$n}->[7] = "$contigs2k";
			    $found++;
			}
			if ($reads2k) {
			    $hash->{$n}->[9] = "$reads2k";
			    $found++;
			}
			$i--;
		    }
		} else {
		    warn "Can't open $notes: $!\n";
		}
	    } else {
		warn "$s_proj has no projectpath or is unavailable\n";
	    }
	    $new_count++;
	    my @array = ( $hash->{$n}->[2], $hash->{$n}->[6], $hash->{$n}->[4], 
			  $hash->{$n}->[1], $hash->{$n}->[11],
			  int($hash->{$n}->[5] + 0.5), 
			  $hash->{$n}->[9], $hash->{$n}->[7], $hash->{$n}->[10]);
	    $table2->add(@array);
	}
    }
    $sth8->finish();

    $dbh->disconnect;

    if ($count == 0){
	print STDOUT "No entries found !\n";
    } else {
	print STDOUT $table2->format();
    }
    warn scalar localtime(), ": Done\n";

}

