#!/usr/local/bin/perl5.6.1 -w

use lib '/nfs/disk100/humpub/modules/PerlModules';

### awol_plate

use strict;
use Symbol ();
use Programs 'gap2caf', 'templates';
use Hum::Tracking 'find_project_directories';
use Getopt::Long 'GetOptions';

{
    my( $major, $cutoff );
    my $fancy = 1;
    my $plain = 0;
    my $write = 0;
    my $help  = 0;
    GetOptions(
        'major=s'   => \$major,
        'cutoff=s'  => \$cutoff,
        'fancy!'    => \$fancy,
        'plain!'    => \$plain,
        'write!'    => \$write,
        'h|help'    => \$help,
        ) or usage();
    usage() if $help;
        
    sub usage {
        exec('perldoc', $0);
    }

    $fancy = 0 if $plain;

    my @projects = @ARGV;
    
    my $proj_dir = find_project_directories(@projects);
    foreach my $proj (@projects) {
        local *STDOUT if $write;
        if ($write) {
            my $report_file = "$proj.awol";
            die "Error: file '$report_file' exists" if -e $report_file;
            open STDOUT, "> $report_file"
                or die "Can't write to '$report_file' : $!";
        }
    
        my $dir = $proj_dir->{$proj}
            or warn "No online path for project '$proj'\n"
            and next;
        #warn "Now in: $dir\n";
        my $pipe = open_template_pipe($proj, $dir);
        
        my $pc = PlateContig->new;
        $pc->majority_fraction($major) if $major;
        $pc->plate_reads_cutoff($cutoff) if $cutoff;
        $pc->process_template_output($pipe, $dir);
        my @suspect = $pc->list_plates_below_cutoff;
        if ($fancy) {
            print "\nReport for project $proj\n", $pc->fancy_report;
            if (@suspect) {
                print "\nSuspect plates in this project:\n";
            } else {
                print "\nNo suspect plates found\n";
            }
        }
        print map "$proj\t$_\n", @suspect;
    }
}

sub open_template_pipe {
    my $project = uc shift;
    my $dir     = shift;
    
    my $fh = Symbol::gensym();
    my $pipe = "cd $dir; gap2caf -project $project -version 0 -bayesian -cutoff 2 2>/dev/null | templates 2>/dev/null |";
    open $fh, $pipe or die "Can't open pipe '$pipe' : $!";
    return $fh;
}

{
    package PlateContig;
    
    sub new {
        my( $pkg ) = @_;
        
        return bless {}, $pkg;
    }
    
    sub process_template_output {
        my( $pc, $pipe, $dir ) = @_;
        
		my( %exclude_plates ) = ();
		my( $i ) = 1;
		while(-e "$dir/readraid_fofn\.$i") {
			next unless open IN,"$dir/readraid_fofn\.$i";
			while (<IN>) {
				chomp;
				my ($well) = $_;
				$well =~ s/\..+//;
				next unless ($well =~ s/[a-z]\d{1,2}$//);
				$exclude_plates{$well} = 1;
			}
			close IN;
			$i++;
		}
        $pc->delete_cutoff_list_cache;
        
        my $read_count = 0;
        while (<$pipe>) {
            #print;
            my( $contig, $suf, $well, @contig_pos ) = (split)[0,1,2,4,5];
            if(!($contig =~ s/^Contig//)) {
                print STDERR "Unexpected contig name format '$contig'\n";
				next;
			}
            if(!($well =~ s/[a-z]\d{1,2}$//)){
                print STDERR "Unexpected well name format '$well'\n";
				next;
			}
			next if(exists $exclude_plates{$well});
            if(!($suf =~ s/.+\.//)) {
                print STDERR "Unexpected suf name format '$suf'\n";
				next;
			}

			$well = "$suf.$well";
            $pc->count_contig_plate($contig, $well);
            $pc->count_plate_contig($well, $contig);
            $pc->contig_length($contig, @contig_pos);
            $read_count++;
        }
        close($pipe) or die "Error from template pipe : exit($?)";
        die "Didn't get any reads" unless $read_count;
        #print STDERR "Got $read_count reads\n";
    }
    
    sub contig_length {
        my( $pc, $contig, @pos ) = @_;
        
        if (@pos) {
            my $length = $pc->{'_contig_length'}{$contig} || 0;
            foreach my $l (@pos) {
                if ($l > $length) {
                    $pc->{'_contig_length'}{$contig} = $l;
                    $length = $l;
                }
            }
        }
        return $pc->{'_contig_length'}{$contig};
    }
    
    sub count_contig_plate {
        my( $pc, $contig, $plate ) = @_;
        
        $pc->{'_contig_plate'}{$contig}{$plate}++;
    }
    
    sub count_plate_contig {
        my( $pc, $plate, $contig ) = @_;
        
        $pc->{'_plate_contig'}{$plate}{$contig}++;
    }
    
    sub majority_fraction {
        my( $pc, $cutoff ) = @_;
        
        if (defined $cutoff) {
            die "Invalid cutoff '$cutoff'" unless $cutoff <= 1;
            $pc->{'_majority_fraction'} = $cutoff;
        }
        $cutoff ||= $pc->{'_majority_fraction'};
        return defined($cutoff) ? $cutoff : 0.6;
    }
     
    sub plate_reads_cutoff {
        my( $pc, $cutoff ) = @_;
        
        if (defined $cutoff) {
            die "Invalid cutoff '$cutoff'" unless $cutoff <= 1;
            $pc->{'_plate_read_cutoff'} = $cutoff;
        }
        $cutoff ||= $pc->{'_plate_read_cutoff'};
        return defined($cutoff) ? $cutoff : 0.3;
    }
   
    sub column_width {
        my( $pc, $width ) = @_;
        
        if (defined $width) {
            die "Invalid column width '$width' (must be at least 3)"
                unless $width =~ /^\d+$/ and $width > 2;
            $pc->{'_column_width'} = $width;
        }
        $width ||= $pc->{'_column_width'};
        return defined($width) ? $width : 4;
    }
    
    sub plate_list {
        my( $pc ) = @_;
        
        return sort keys %{$pc->{'_plate_contig'}};
    }
    
    sub plate_contig_list {
        my( $pc, $plate ) = @_;
        
        return sort keys %{$pc->{'_plate_contig'}{$plate}};
    }
    
    sub contig_list {
        my( $pc ) = @_;
        
        return sort {$a <=> $b} keys %{$pc->{'_contig_plate'}};
    }
    
    sub contig_list_by_read_count {
        my( $pc ) = @_;
        
        return map $_->[0], $pc->contig_read_count_list;
    }
    
    sub contig_read_count_list {
        my( $pc ) = @_;
        
        my( @sort );
        foreach my $contig ($pc->contig_list) {
            push(@sort, [$contig, $pc->contig_read_count($contig)]);
        }
        return sort {$b->[1] <=> $a->[1] or $b->[0] cmp $a->[0]} @sort;
    }
    
    sub contig_read_count_cutoff_list {
        my( $pc ) = @_;
        
        my( $cutoff_list );
        unless ($cutoff_list = $pc->{'_contig_read_count_cutoff_list'}) {
            $cutoff_list = [];
            my $cutoff = $pc->majority_fraction;
            my $total = $pc->total_read_count
                or die "Total contig count returned 0";
            my $num_contigs_seen = 0;
            foreach my $cc ($pc->contig_read_count_list) {
                push(@$cutoff_list, $cc->[0]);
                $num_contigs_seen += $cc->[1];

                # Stop adding plates to the list when we reach
                # the majority_fraction
                last if ($num_contigs_seen / $total) > $cutoff;
            }
            $pc->{'_contig_read_count_cutoff_list'} = $cutoff_list;
        }
        return @$cutoff_list;
    }
    
    sub delete_cutoff_list_cache {
        my( $pc ) = @_;
        
        $pc->{'_contig_read_count_cutoff_list'} = undef;
    }
    
    sub proportion_plate_within_cutoff {
        my( $pc, $plate, $report_total ) = @_;
        
        my %is_majority_contig = map {$_, 1} $pc->contig_read_count_cutoff_list;
        my($within, $without) = (0,0);
        foreach my $contig ($pc->plate_contig_list($plate)) {
            my $count = $pc->plate_contig_count($plate, $contig);
            if ($is_majority_contig{$contig}) {
                $within += $count;
            } else {
                $without += $count;
            }
        }
        my $total = $within + $without
            or die "No reads found for plate '$plate'";
        return ($within / $total,$total) if $report_total;
        return ($within / $total);
    }
    
    sub list_plates_below_cutoff {
        my( $pc ) = @_;
        
        my $cutoff = $pc->plate_reads_cutoff;
        my( @below_cutoff );
        foreach my $plate ($pc->plate_list) {
			my ($pcv,$total) = $pc->proportion_plate_within_cutoff($plate,1);
            push(@below_cutoff, ((split /\./,$plate)[1].".".(split /\./,$plate)[0]))
                if $pcv < $cutoff && $total > 10;
        }
        return @below_cutoff;
    }
    
    sub contig_read_count {
        my( $pc, $contig ) = @_;
        
        my $total = 0;
        foreach my $plate ($pc->plate_list) {
            $total += $pc->plate_contig_count($plate, $contig);
        }
        return $total;
    }
    
    sub total_read_count {
        my( $pc ) = @_;
        
        my $total = 0;
        my @contigs = $pc->contig_list;
        foreach my $plate ($pc->plate_list) {
            foreach my $contig (@contigs) {
                $total += $pc->plate_contig_count($plate, $contig);
            }
        }
        return $total;
    }
    
    sub contig_plate_count {
        my( $pc, $contig, $plate ) = @_;
        
        if (exists $pc->{'_contig_plate'}{$contig}{$plate}) {
            return $pc->{'_contig_plate'}{$contig}{$plate};
        } else {
            return 0;
        }
    }
    
    sub plate_contig_count {
        my( $pc, $plate, $contig ) = @_;
        
        if (exists $pc->{'_plate_contig'}{$plate}{$contig}) {
            return $pc->{'_plate_contig'}{$plate}{$contig};
        } else {
            return 0;
        }
    }
    
    sub fancy_report {
        my( $pc ) = @_;
        
        my($header, @plates) = $pc->fancy_plate_header;
        
        my $cw = $pc->column_width;
        
        my $table = '';
        {
            my @contigs = $pc->contig_list;
            my $contigs_pattern = "%5s %6d ". ("\%${cw}s" x @plates) ."\n";
            foreach my $contig (@contigs) {
                my $contig_length = $pc->contig_length($contig);
                my( @counts );
                foreach my $plate (@plates) {
                    push(@counts, $pc->contig_plate_count($contig, $plate) || '.');
                }
                $table .= sprintf($contigs_pattern, $contig, $contig_length, @counts);
            }
        }
        
        my $percent_pattern = '   %% in major' . ("\%${cw}.0f" x @plates) ."\n";
        my $percent_in_majority_line = sprintf($percent_pattern,
            map $pc->proportion_plate_within_cutoff($_,0) * 100, @plates);
        
        my $cutoff_line = sprintf("Majority fraction = %2.0f%%  Plate cutoff = %2.0f%%\n",
            $pc->majority_fraction  * 100, 
            $pc->plate_reads_cutoff * 100);
        
        #my $total_read_line = "Total read count = ". $pc->total_read_count ."\n";
        return join("\n", $header, $table, $percent_in_majority_line, $cutoff_line);
    }
    
    sub fancy_plate_header {
        my( $pc, $left_offset ) = @_;
        
        my $cw = $pc->column_width;
        
        my $left_padding = ' ' x ($left_offset || 13);
        my @plates = $pc->plate_list;
        #my @plates = (qw{bA11M8-1 bA12M8-2 bA21M8-3 bA21M8-34});
        my( @sort );
        foreach my $t (@plates) {
            my( $proj, $num ) = $t =~ /(.+?)\W*(\d+)$/
                or die "Unexpected plate name format '$t'";
            push(@sort, [$t, $proj, $num]);
        }
        @sort = sort {$a->[1] cmp $b->[1] or $a->[2] <=> $b->[2]} @sort;
        @plates = map $_->[0], @sort;
        my @plate_line;
        my $dash_seg            = '-' x $cw;
        my $gap_dash_seg        = '  '. ('-' x ($cw - 2));
        my $number_line         = '';
        my $dash_line           = '';
        for (my $i = 0; $i < @sort; $i++) {
            my($t, $proj, $num) = @{$sort[$i]};
            $number_line .= sprintf("\%${cw}s", $num);
            if ($i == 0 or $sort[$i-1][1] ne $proj) {
                #print STDERR "New proj $proj\n";
                my $header_offset = $i * $cw;
                for (my $j = 0;; $j++) {
                    my $line = $plate_line[$j] || '';
                    my $padding = $header_offset - length($line);
                    #print STDERR "Padding = $padding\n";
                    if ($padding >= 0) {
                        # We need at least 1 space, unless it is the first line
                        next if $padding == 0 and $i != 0;
                        $line .= (' ' x $padding) . " $proj";
                        $plate_line[$j] = $line;
                        last;
                    }
                }
                $dash_line .= $gap_dash_seg;
            } else {
                $dash_line .= $dash_seg;
            }
        }
        @plate_line = map " $_", reverse @plate_line;
        my $head = "\n" . join('', map "$left_padding$_\n", (@plate_line, $dash_line, $number_line));
        return $head, @plates;
    }
}

__END__

=head1 NAME - awol_plate

=head1 SYNOPSIS

    # Investigate projects bA21M8 and bA29C8
    # for awol plates:
    awol_plate -write bA21M8 bA29C8

=head1 OPTIONS

=over4

=item -cutoff <FRACTION>

Specify the fraction of reads from a plate which
have been assembled into the contigs containing
the majority of reads below which a plate
considered suspect.  Default is B<0.3>.

=item -major <FRACTION>

Specify the fraction of the total number of reads
used to decide which contigs constitute the
majority.  Default is B<0.6>.

=item -fancy

Print a full human-readable report, plus the list
output by the C<-plain> option below.

=item -plain

Just output a tab-delimited list of PROJECT and
PLATE names.

=item -write

Write the report for each project to a file named
C<PROJECT.awol>.

=back

=head1 DESCRIPTION

B<awol_plate> is designed to spot plates which do not
belong in projects.  At its default settings, it
will report plates where less than 30% of the
reads from the plate contribute to the contigs
which contain 60% of the total number of reads in
the project.

=head2 The fancy format

Plates are listed along the header of the report
table.

The first column is the number of the contig from
the gap database, and the second is the length of
the contig.

The body of the table show how many reads in each
contig came from which plate.

The footer of the table shows the percentage of
reads from each plate which contribute to contigs
containing the majority of reads, and shows the
cutoffs used.

=head1 AUTHOR

James Gilbert B<email> jgrg@sanger.ac.uk

=head1 HISTORY

 # $Log: awol_plate.pl,v $
 # Revision 1.2  2003/03/19  17:21:17  rmd
 # Now uses perl5.6.1
 #
 # 
 # Revision 1.1  2000/08/18  14:34:43  jcm
 # Initial revision

=cut
