# -*-Perl-*- 
#
# $Id: GNUtils.pm,v 1.114 2006/10/25 15:54:19 js10 Exp $
#
#

package GNUtils;
use WrapDBI;
use Caftools;
use GD;
use Exporter ();
use Cwd;

use Data::Dumper;

@ISA = qw(Exporter);
@EXPORT = qw(find_problem_regions);

use strict;

sub new {
  my($type,$caf_obj,$checkingmode) = @_;

  my $dir = getcwd();
  $dir =~ s/.*\///;
  my $self = {
      'CAF'     => $caf_obj,
      'CMODE'   => $checkingmode, # cmode is set by auto checking.
      'PROJECT' => $dir          
      };
  
  bless $self, $type;
  
  return $self; 
}

sub find_problem_regions {
  my ($self,$searches) = @_;
  
  $self->parse_args($searches);
  
  my(%raid,%wraid,%addi,%sil,%pcr,@problems,%seen,%sil_origin,%consensus,%wgs,%parent_sil);
  
  #
  # Get hash list of suffixes from oracle.
  # also fill the sil parent and pcr hashes.
  my $suffix = $self->init_suffix_types(\%sil_origin,\%addi,\%pcr,\%parent_sil);
  

  #
  # Obtain hash lists for any 'special' regions ...
  #
  foreach my $read ($self->{CAF}->readsList) {
    my $read_info = $self->{CAF}->getEntry($read);
    my $read_data = \$read_info->{Sequence};
    my $read_qual = \$read_info->{BaseQuality};
    my %qual_hash;
    map { $qual_hash{$_}++ } (split" ",$$read_qual);
    my $q_type_count += keys %qual_hash;

    # find consensus style reads.
    if(($$read_data =~ /Clone unknown/) &&
       ($$read_data =~ /Sequencing_vector "unknown"/) &&
       ($$read_data !~ /wwRaided/)&& $q_type_count == 1) {
        $consensus{$read}++;
    }
    # Read raid reads.
    if($$read_data =~ /Stolen/) {
        if($self->{RR}) { $raid{$read}++ if (($$read_data !~ /wwRaided/) && ($read !~ /c\d+\.Contig.*$/) && ($read !~ /^contig_\d+/)); }
        if($self->{WW}) { $wraid{$read}++ if ($$read_data =~ /wwRaided/); }        
        if($self->{WGS}) { $wgs{$read}++ if (($read =~ /c\d+\.Contig.*$/) || ($read =~ /^contig_\d+/)); }
    }
    if($$read_data =~ /FWGS/) { $wgs{$read}++; } # For new tag type but some are still STOL..so include the above.
    if($self->{PCR}) { $pcr{$read}++ if ($read =~ /\.w/); }
  }
  
  #
  # Iterate through the contig lists finding regions
  # and checking accordingly
  #
  foreach my $contig ($self->{CAF}->contigsList()) {
      my $contig_entry = $self->{CAF}->getEntry($contig);
      my ($assembled_from,$strand,$last,$eoc) = $self->get_assembled_from(\$contig_entry->{Sequence});  
      # Skip single read contigs.
      if (scalar(@{$assembled_from})==2) { next; }
      next unless (@{$assembled_from});
      my $sv = 0;
      my %region=();
      
      #foreach contig in the database...
      while(@{$assembled_from}){

          # fetch the next region. (defined as a run between one reads end and another reads start).
          $self->fetch_region(\$sv,\@{$assembled_from},\%region);                 
          &remove_consensus_reads(\%consensus,\%region);
          next if (!scalar %region);
          my $ac += keys %region;          
          $self->region_checks(\%region,$ac,$suffix,\%addi,\%pcr,
                               \%raid,\%wraid,$strand,$sv,$contig,\%seen,
                               $$assembled_from[0]->[1],\@problems,$last,$eoc,\%sil_origin,\%wgs,\%parent_sil);
          
      }
      # Single functions that can run on the entire contig.
      #
      # Screen for low quality bases.
      if($self->{QUAL}) { $self->screen_low_quality($contig,\@problems,\$contig_entry->{DNA},\$contig_entry->{BaseQuality}); }

      # Single base pair runs. 10 or more of any single base....
      if($self->{SBPR}) { $self->single_base_pair_runs(\@problems,\$contig_entry->{DNA},$contig); }
      
      # TA run repeats.
      if($self->{TAREP}) { $self->ta_repeat_screen(\@problems,\$contig_entry->{DNA},$contig); }

      # low qual wgs regions.
      if($self->{WGSS}) { 	  
	  my $wgslowqual = $self->screen_low_quality_wgs_data($contig,\$contig_entry->{Sequence}); 
	  foreach my $prob_ref (@problems) {
	      next unless ($prob_ref->[3]=~/Uni\-directional|MouseWGS\sdata\swith\s1|PCR|SIL|Single\stemplate|ReadRaid\sonly/);
	      for(my $i = $prob_ref->[1]+1; $i<= $prob_ref->[2]; $i++) {
                  if(exists $$wgslowqual{$i}) {
		      push @problems, [$$wgslowqual{$i}->[0],
				       $$wgslowqual{$i}->[1],
				       $$wgslowqual{$i}->[1],
				       $$wgslowqual{$i}->[2],
				       $i
				       ];				       
		  }
	      }	      
	  }
      }
  }
  
  # Sort problems list by position. 

  # added 1/9/06. HACK. 
  # for SIL -> retrans (Single:SIL) need to reparse regions
  # to see if SIL regions are next to each other, in which case join them
  #
  # 
  #  single template ------- 
  #              SIL ------- ------- SIL  (JOIN)
  $self->sort_list(\@problems);
  
  return \@problems;
  
}

sub remove_consensus_reads {
    my($cons_ref,$read_ref) = @_;    
    foreach my $conseq ( keys %$cons_ref )  {
        if(exists $$read_ref{$conseq}) {
            delete $$read_ref{$conseq};
        }
    }
}

#
# add leading 0 to all well numbers.
#
sub fix_read {
  my($self,$read_ref) = @_;  

	$$read_ref =~ s/([a-p])(\d)\./sprintf("%s%02d.", $1, $2)/e;
	$$read_ref =~ s/^([^.]+[a-p])(\d)$/sprintf("%s%02d", $1, $2)/e;
}

sub region_checks {
    my ($self,$rlist,$count,$chem,$addi,$pcr,$raid,$wraid,$strand,$start,$contig,$seen,$end,$prob,$last,$eoc,$parent,$wgs,$parent_sil) = @_;
    
    # We supply a hash containing reads covering a region which
    # is best described as ...
    # 'an area between one reads end point and anothers start'
    # and perform subroutines on those reads to find the problem areas.
    #
    # N.B $count is the count of the number of reads in this region.

    # These follow the same format.
    # 1) see if we have a problem region
    # 2) check if start for problem is already set. if it is :next else set start for problem.
    # 3) keep :next'ing while we have a problem region and start is defined.
    # 4) no longer in region, store problem region (start - current pos) and undefine start.

    # SEE ST BELOW FOR EXAMPLE.


    if($self->{ST}) {
        # Single Template.
        my %seen=();
        foreach my $name (keys %$rlist) {
            # cope with leading zeros etc...
            $self->fix_read(\$name);
            my ($tp,$sf) = split/\./,$name;
            $seen{$tp}++;      
        }       
        $self->{current_t} = keys %seen;
        # in problem region, start not set.....
        if(!exists $self->{stset} && $self->{current_t} == 1) {
            # print if end of contig.
            if($start == $last){
                $self->gap_format($contig,$start,"Single template",$eoc,$prob); 
                delete $self->{stset};
            }else{
                # problem region....set start....
                $self->{stset} = $start;
            }              
            # start set but out of problem region.... so report..
        }elsif(exists $self->{stset} && $self->{current_t} != 1) {
            $self->gap_format($contig,$self->{stset},"Single template",$start,$prob);
            delete $self->{stset};
            
            # start set, still in problem region but at end of contig.... so report...
        }elsif(exists $self->{stset} && $self->{current_t} == 1 && $start == $last) {
            $self->gap_format($contig,$self->{stset},"Single template",$eoc,$prob);
            delete $self->{stset};
            # do nothing ... keep iterating ....
        }else{}
    }
    
    if($self->{'2R'}) {
        # + Only two reads
        if(!exists $self->{tworeadset} && $count==2) {
            if($start == $last){
                $self->gap_format($contig,$start,"Only two reads",$eoc,$prob);
                delete $self->{tworeadset};
            }else{
                $self->{tworeadset} = $start;
            }
        }elsif(exists $self->{tworeadset} && $count != 2) {
            $self->gap_format($contig,$self->{tworeadset},"Only two reads",$start,$prob);
            delete $self->{tworeadset};
        }elsif(exists $self->{tworeadset} && $count==2 && $start == $last) {
            $self->gap_format($contig,$self->{tworeadset},"Only two reads",$eoc,$prob);
            delete $self->{tworeadset};
        }else{}
    }
    
    if($self->{'3R'}) {
        # + Only three reads
        if(!exists $self->{threereadset} && $count==3) {
            if($start == $last) {
                $self->gap_format($contig,$start,"Only three reads",$eoc,$prob);
                delete $self->{threereadset};
            }else{
                $self->{threereadset} = $start;
            }	
        }elsif(exists $self->{threereadset} && $count!=3) {
            $self->gap_format($contig,$self->{threereadset},"Only three reads",$start,$prob);
            delete $self->{threereadset};
        }elsif(exists $self->{threereadset} && $count==3 && $start == $last) {
            $self->gap_format($contig,$self->{threereadset},"Only three reads",$eoc,$prob);
            delete $self->{threereadset};
        }else{}
    }
    
    if($self->{ADDI}) {
        # + Additive regions only
        my $acount=0;
        my $count_t = $count;


        #
        # This occurs throught the searches. just removes the WGS (whole genome shotgun) reads
        # from the region asn checking didnt want them included.
        #
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
            if(exists $$addi{$name}) {$acount++;}
        }
        # Add this to all others to fix bug where a region containing wgs only (which 
        # is ignored) will show up as all the other regions. read count = 0, problem count = 0
        # It will always succeed when there is nothing there.
        if($count_t == 0) {$count_t=1; }

        if(!exists $self->{addiset} && $acount == $count_t) {
            if($start == $last){
                $self->gap_format($contig,$start,"Additive only",$eoc,$prob);	
                delete $self->{addiset};
            }else{
                $self->{addiset} = $start;
            }	
        }elsif(exists $self->{addiset} && $acount != $count_t) {
            $self->gap_format($contig,$self->{addiset},"Additive only",$start,$prob);
            delete $self->{addiset};
        }elsif(exists $self->{addiset} && $acount == $count_t && $start == $last && $count_t > 0) {
            $self->gap_format($contig,$self->{addiset},"Additive only",$eoc,$prob);
            delete $self->{addiset};
        }else{}
    }

    if($self->{UADDI}) {
        # + Uni-directional ADDI only regions, either forward or reverse.
        my $adrev_strandc  = 0;
        my $adfwd_strandc  = 0;
        my $acount=0;
        my $pcount=0;
        my $udtype = undef;
        my $count_t = $count;
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
		if($count > 1) {$count_t -= 1;}
                next;
            }
			my ($suf) = $name =~ /\.(...)/;
            if(exists $$addi{$name}) {$acount++;}
            if($suf && exists $$chem{primer}->{$suf}) {$pcount++;}
            if(exists $$strand{$name}) {$adrev_strandc++;}
            if(!exists $$strand{$name}) {$adfwd_strandc++;}
        }      
        # Need to include primers within the region as these will need to 
        # be tagged. Also need to work out when the UADDI region (possibly)
        # moves to primers only and format the previous region, if defined.
        if($count_t == 0) {$count_t=1; }
        if($pcount >= $count_t && exists $self->{unidirection_addi}) {
            # UADDI Region already defined...
            $self->gap_format($contig,$self->{unidirection_addi},$self->{last_addi_strand_type},$start,$prob);
            delete $self->{unidirection_addi};
            delete $self->{last_addi_strand_type};
        }else{
            # add primer count to additive count to help the next bit work.
            $acount += $pcount;
            
            if(($adfwd_strandc == $acount)&&($acount==$count_t)) { $udtype = 'Uni-directional Additive only (Forward strand)'; }
            if(($adrev_strandc == $acount)&&($acount==$count_t)) { $udtype = 'Uni-directional Additive only (Reverse strand)'; }      
            
            if(!exists $self->{unidirection_addi} && (defined $udtype)) {
                if($start == $last) {
                    $self->gap_format($contig,$start,$udtype,$eoc,$prob);
                    delete $self->{unidirection_addi};
                    delete $self->{last_addi_strand_type};
                }else{
                    $self->{unidirection_addi} = $start;                
                    $self->{last_addi_strand_type} = $udtype;
                }	
            }elsif(exists $self->{unidirection_addi} && (!defined $udtype)) {
                $self->gap_format($contig,$self->{unidirection_addi},$self->{last_addi_strand_type},$start,$prob);
                delete $self->{unidirection_addi};
                delete $self->{last_addi_strand_type};
            }elsif(exists $self->{unidirection_addi} && (defined $udtype && $udtype ne $self->{last_addi_strand_type})){
                $self->gap_format($contig,$self->{unidirection_addi},$self->{last_addi_strand_type},$start,$prob);
                $self->{unidirection_addi} = $start;
                $self->{last_addi_strand_type} = $udtype;
            }elsif(exists $self->{unidirection_addi} && (defined $udtype) && ($start == $last)) {
                $self->gap_format($contig,$self->{unidirection_addi},$self->{last_addi_strand_type},$eoc,$prob);
                delete $self->{unidirection_addi};
                delete $self->{last_addi_strand_type};
            }else{}
        } 
    }
    

    #
    # Pain and suffering here Pt.1 !! A multitude of checks required some 
    # tricky SQL. mostly to do with finding the parent (however many levels) of
    # each read and checking whether it's over the region...with its sibling...without etc...
    if($self->{SIL}) {
        # + SIL regions only
        my $scount = 0;  
        my $count_t = $count;
        my $til_count = 0;
        my $retrans_count = 0;
        my $type;
        my (%parent_count,%seen,%types);
        if($count_t == 0) {$count_t=1; }
        # ignore WGS reads, remove from count
        foreach my $name (keys %$rlist) {	        
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
        }
        # Need to count. one for number of sils (+corresponding parents as they are sil)
        # ...and another that counts the number of original parents over this region.
        foreach my $name (keys %$rlist) {	            
            my $target;
            if($name =~ /^([J|T])(\d+)(\-)?(\S+)?[a-p]\d{1,2}\..*$/) {                
                if($1 eq 'T') { $til_count++; } # count the tils for checking later.
                $target = $2;
                $scount++ if (exists $parent->{$target}); # is a sil read, is in the hash. If RFSIL, wont be found.
                my $parent_name = $parent->{$target} || undef;            
                
                if(defined $parent_name) {
                    $parent_count{$parent_name}++;
                    # find out if parent(s) is over the region and return how many.
                    # count these as sil.
                    if(!$seen{$parent_name}){   
                        $scount += &parent_over_region($rlist,$parent_name);
                        # ...have parent(s) then count...
                        if($scount>0) {
                            $seen{$parent_name}++;
                        }
                    }                    
                }

            # Retransformation bit.
            }elsif($name =~ /^([R]\d+(\w)?\-?\S+?[a-p]\d{1,2})\..*$/) {
                $retrans_count++; # count the retrans for checking later.
                $target = $1;
                $scount++; # is a sil read.
                my $parent_name = $parent->{$target} || undef; # get my parent.
                
                if(defined $parent_name) {
                    $parent_count{$parent_name}++; # count unique parents.

                    # find out if parent(s) is over the region and return how many.
                    # count these as sil.
                    if(!$seen{$parent_name}){
                        # Find parents, return number found.
                        $scount += &parent_over_region($rlist,$parent_name);
                        # ...have parent(s) then count as looked for...
                        if($scount>0) {
                            $seen{$parent_name}++;
                        }
                    }
                }                              
                # To save multiple hacks below i'll set this to SIL:Single
                # and split into the 2 regions within 'sub gap_format'.
                if($scount == $count_t && $retrans_count){
                    $type = "SIL:Single";
                }
            }           
        }
               
        # NB:to fix: these are the same. just use one hash.
        my $total_parents += keys %parent_count;              
        my $total_parents_over_region += keys %seen;

        # Need to work out type if not defined above.       
        if(!defined $type) {  $type = 'SIL only'; } 

        # if TILs only...
        if($til_count == $count_t) { $type = 'TIL only'; }
       
 
        if(!exists $self->{silset} && $scount == $count_t && $total_parents < 2) {
            if($start == $last){                
                $self->gap_format($contig,$start,$type,$eoc,$prob);
                delete $self->{silset};
            }else{
                $self->{silset} = $start;
                $self->{last_sil_type} = $type;
            }      
            # find those regions that are sil only BUT from single parents. multiples are double-stranded.
        }elsif(exists $self->{silset} && $scount == $count_t && ($total_parents > 1 || $type ne $self->{last_sil_type})){
            $self->gap_format($contig,$self->{silset},$self->{last_sil_type},$start,$prob);                      
            delete $self->{silset}; 
            if($type ne $self->{last_sil_type}) {
                $self->{last_sil_type} = $type; # set to new type. 
                $self->{silset} = $start;       # NB: reset silset to current start as we're still in a SIL region.
            }else{
                delete $self->{last_sil_type};
            }

        }elsif(exists $self->{silset} && ($scount != $count_t || $count_t == $total_parents_over_region)){
            $self->gap_format($contig,$self->{silset},$self->{last_sil_type},$start,$prob);
            delete $self->{silset};
            delete $self->{last_sil_type};            

        }elsif(exists $self->{silset} && $scount == $count_t && $start == $last) {            
            $self->gap_format($contig,$self->{silset},$self->{last_sil_type},$eoc,$prob);
            delete $self->{silset};
            delete $self->{last_sil_type}; 
        }else{}

        
    }


    if($self->{PCR}) {
        # + PCR regions only
        my $pcount = 0;    
        my $count_t = $count;
        my %parent_count;
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
            # find PCRs, skip clone walks (starting W)
            if(exists $$pcr{$name} && $name !~ /^W/) {$pcount++;}            
            
            # Pick up those that are PCR for SIL.
            if($name =~ /^[J|T|S]?(\d+(\w))?(\-)?(\S+)?[a-p]\d{1,2}..*$/) {
                if(exists $$pcr{$1}) {$pcount++;}  
            }
        }        
        if($count_t == 0) {$count_t=1; }
        if(!exists $self->{pcrset} && $pcount == $count_t) {
            if($start == $last) {
                $self->gap_format($contig,$start,"PCR only",$eoc,$prob);
                delete $self->{pcrset};
            }else{
                $self->{pcrset} = $start;
            }
        }elsif(exists $self->{pcrset} && $pcount != $count_t) {
            $self->gap_format($contig,$self->{pcrset},"PCR only",$start,$prob);
            delete $self->{pcrset};
        }elsif(exists $self->{pcrset} && $pcount == $count_t && $start == $last) {
            $self->gap_format($contig,$self->{pcrset},"PCR only",$eoc,$prob);
            delete $self->{pcrset};
        }else{}
    }
    
    
    # Special case. For Readraid, WWraid and WGS we need to find all regions containing
    # 'raid + <readnum>' or less, but we need to write a problem line everytime
    # the regions changes. i.e Readraid + 1 read *STOP* Readraid only *STOP* ...etc..
    # Also, allow a user defined number of templates (formerly project reads).
    
    if($self->{RR}) {
        # + Read raid covered regions
        my $rrcount = 0;
        my $nrcount = 0;   # read count;
        my $templates = 0; # template count.
        my $tp;
        my %seen=();
        my $count_t = $count;
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
            if(exists $$raid{$name}) { 
                $rrcount++; 
            }else{
                $self->fix_read(\$name);
                my ($tp) = (split/\./,($name))[0];
                ++$seen{$tp};
                ++$nrcount;
            }
        }    
        if($count_t == 0) {$count_t=1; }
        $templates = keys %seen;       
        if(!exists $self->{rrset} && $templates <= $self->{R} && $count_t == ($rrcount+$nrcount) && $count_t != $nrcount) {
            if($start == $last) {
                if($nrcount) { $tp = "ReadRaid with $nrcount template(s)"; }else{ $tp = "ReadRaid only";}
                $self->gap_format($contig,$start,$tp,$eoc,$prob);
                delete $self->{rrset};
            }else{
                $self->{rrset} = $start;      
                $self->{nrset} = $templates;
            }	
        }elsif(exists $self->{rrset} && ($templates > $self->{R} || $count_t == $nrcount || $count_t != ($rrcount+$nrcount))) {
            if(exists $self->{nrset} && $self->{nrset} > 0) { $tp = "ReadRaid with $self->{nrset} template(s)"; }else{ $tp = "ReadRaid only";}
            $self->gap_format($contig,$self->{rrset},$tp,$start,$prob);
            delete $self->{rrset};
            delete $self->{nrset};
        }elsif(exists $self->{rrset} && $templates <= $self->{R} && $count_t == ($rrcount+$nrcount) && $count_t != $nrcount) {
            if($start == $last) {
                if(exists $self->{nrset} && $self->{nrset} > 0) { $tp = "ReadRaid with $self->{nrset} template(s)"; }else{ $tp = "ReadRaid only";}
                $self->gap_format($contig,$self->{rrset},$tp,$eoc,$prob);
                delete $self->{rrset};
                delete $self->{nrset};
            }elsif($templates != $self->{nrset}){
                if(exists $self->{nrset} && $self->{nrset} > 0) { $tp = "ReadRaid with $self->{nrset} template(s)"; }else{ $tp = "ReadRaid only";}
                $self->gap_format($contig,$self->{rrset},$tp,$start,$prob);
                $self->{rrset} = $start;  
                $self->{nrset} = $templates;
            }else{}
        }else{}
    }
    
    if($self->{WW}) {
        # + WWraid covered regions
        my $wcount = 0;
        my $nwcount = 0;   # read count.
        my $templates = 0; # template count.
        my $tp;
        my %seen;
        my $count_t = $count;
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
            if(exists $$wraid{$name}) {
                $wcount++;            
            }else{
                $self->fix_read(\$name);
                my ($tp) = (split/\./,($name))[0];
                ++$seen{$tp};
                ++$nwcount;
            }
        }     
        if($count_t == 0) {$count_t=1; }
        $templates = keys %seen;
        if(!exists $self->{wrset} && ($templates <= $self->{W}) && ($count_t == ($wcount+$nwcount)) && ($count_t != $nwcount)) {
            if($start == $last) {
                if($nwcount) { $tp = "WWRaid with $nwcount template(s)"; }else{ $tp = "WWRaid only";}
                $self->gap_format($contig,$start,$tp,$eoc,$prob);
                delete $self->{wrset};
            }else{
                $self->{wrset} = $start;      
                $self->{nwset} = $templates;
            }	
        }elsif(exists $self->{wrset} && ($templates > $self->{W} || $count_t == $nwcount || $count_t != ($wcount+$nwcount))) {
            if(exists $self->{nwset} && $self->{nwset} > 0) { $tp = "WWRaid with $self->{nwset} template(s)"; }else{ $tp = "WWRaid only";}
            $self->gap_format($contig,$self->{wrset},$tp,$start,$prob);
            delete $self->{wrset};
            delete $self->{nwset};
        }elsif(exists $self->{wrset} && $templates <= $self->{W} && $count_t == ($wcount+$nwcount) && $count_t != $nwcount) {
            if($start == $last) {
                if(exists $self->{nwset} && $self->{nwset} > 0) { $tp = "WWRaid with $self->{nwset} template(s)"; }else{ $tp = "WWRaid only";}
                $self->gap_format($contig,$self->{wrset},$tp,$eoc,$prob);
                delete $self->{wrset};
                delete $self->{nwset};
            }elsif($templates != $self->{nwset}){
                if(exists $self->{nwset} && $self->{nwset} > 0) { $tp = "WWRaid with $self->{nwset} template(s)"; }else{ $tp = "WWRaid only";}
                $self->gap_format($contig,$self->{wrset},$tp,$start,$prob);
                $self->{wrset} = $start;  
                $self->{nwset} = $templates;
            }else{}
        }else{}
    }
    
    if($self->{WGS}) {
        # + Mouse WGS covered regions
        my $wgscount = 0;
        my $nwgscount = 0; # read count.
        my $templates = 0; # template count.
        my $tp;
        my %seen;
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                $wgscount++;
            }else{
                $self->fix_read(\$name);
                my ($tp) = (split/\./,($name))[0];
                ++$seen{$tp};
                ++$nwgscount;
            }
        }  
        $templates = keys %seen;
        if(!exists $self->{wgsset} && ($templates <= $self->{M}) && ($count == ($wgscount+$nwgscount)) && ($count != $nwgscount)) {
            if($start == $last) {
                if($nwgscount) { $tp = "MouseWGS data with $templates template(s)"; }else{ $tp = "MouseWGS only";}
                $self->gap_format($contig,$start,$tp,$eoc,$prob);
                delete $self->{wgsset};
            }else{
                $self->{wgsset} = $start;      
                $self->{nwgsset} = $templates;
            }	
        }elsif(exists $self->{wgsset} && ($templates > $self->{M} || $count == $nwgscount || $count != ($wgscount+$nwgscount))) {
            if(exists $self->{nwgsset} && $self->{nwgsset} > 0) { $tp = "MouseWGS data with $self->{nwgsset} template(s)"; }else{ $tp = "MouseWGS only";}
            $self->gap_format($contig,$self->{wgsset},$tp,$start,$prob);
            delete $self->{wgsset};
            delete $self->{nwgsset};
        }elsif(exists $self->{wgsset} && $templates <= $self->{M} && $count == ($wgscount+$nwgscount) && $count != $nwgscount) {
            if($start == $last) {
                if(exists $self->{nwgsset} && $self->{nwgsset} > 0) { $tp = "MouseWGS data with $self->{nwgsset} template(s)"; }else{ $tp = "MouseWGS only";}
                $self->gap_format($contig,$self->{wgsset},$tp,$eoc,$prob);
                delete $self->{wgsset};
                delete $self->{nwgsset};
            }elsif($templates != $self->{nwgsset}){
                if(exists $self->{nwgsset} && $self->{nwgsset} > 0) { $tp = "MouseWGS data with $self->{nwgsset} template(s)"; }else{ $tp = "MouseWGS only";}
                $self->gap_format($contig,$self->{wgsset},$tp,$start,$prob);
                $self->{wgsset} = $start;  
                $self->{nwgsset} = $templates;
            }else{}
        }else{}
    }


    if($self->{FST}) {
        # + Forward strand (+ve) and terminator only
        my $fst = 0;
        my $strandc = 0;
        my $count_t = $count;
        my $wgs_set = undef;
        if($self->{CMODE}) {
            foreach my $name (keys %$rlist) {
                if(exists $$wgs{$name}) {
                    $wgs_set++;
                }
            }
        } 
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
            # Get the suffix.
			my ($suf) = $name =~ /\.(...)/;
            if($self->{CMODE} && $wgs_set) {
                #don't incrememnt.
            }else{
                if(!exists $$strand{$name}) {$strandc++;}
                if($suf && exists $$chem{term}->{$suf} || exists $$chem{bdt}->{$suf}) {$fst++;}        
            }
        }
        if($count_t == 0) {$count_t=1; }
        if(!exists $self->{fsto} && $fst == $count_t && $strandc == $count_t) {
            if($start == $last) {
                $self->gap_format($contig,$start,"Forward strand and terminator only",$eoc,$prob);
                delete $self->{fsto};
            }else{
                $self->{fsto}= $start;      
            }
        }elsif(exists $self->{fsto} && ($fst != $count_t || $strandc != $count_t)) {
            $self->gap_format($contig,$self->{fsto},"Forward strand and terminator only",$start,$prob);
            delete $self->{fsto};
        }elsif(exists $self->{fsto} && $fst == $count_t && $strandc == $count_t && $start == $last) {
            $self->gap_format($contig,$self->{fsto},"Forward strand and terminator only",$eoc,$prob);
            delete $self->{fsto};
        }else{}    
    }
    
    if($self->{RST}) {
        # + Reverse strand (-ve) and terminator only
        my $rst = 0;
        my $strandc = 0;
        my $count_t = $count;
        my $wgs_set = undef;
        if($self->{CMODE}) {
            foreach my $name (keys %$rlist) {
                if(exists $$wgs{$name}) {
                    $wgs_set++;
                }
            }
        }
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
            # Get the suffix.
			my ($suf) = $name =~ /\.(...)/;
            if($self->{CMODE} && $wgs_set) {
                #don't increment.
            }else{
                if(exists $$strand{$name}) {$strandc++;}
                if($suf && exists $$chem{term}->{$suf} || exists $$chem{bdt}->{$suf}) {$rst++;}        
            }
        }
        if($count_t == 0) {$count_t=1; }
        if(!exists $self->{rsto} && $rst == $count_t && $strandc == $count_t) {
            if($start == $last) {
                $self->gap_format($contig,$start,"Reverse strand and terminator only",$eoc,$prob);
                delete $self->{rsto};
            }else{
                $self->{rsto} = $start;
            }
        }elsif(exists $self->{rsto} && ($rst != $count_t || $strandc != $count_t)) {
            $self->gap_format($contig,$self->{rsto},"Reverse strand and terminator only",$start,$prob);
            delete $self->{rsto};
        }elsif(exists $self->{rsto} && $rst == $count_t && $strandc == $count_t && $start == $last) {
            $self->gap_format($contig,$self->{rsto},"Reverse strand and terminator only",$eoc,$prob);
            delete $self->{rsto};
        }else{}
    }
    
    if($self->{FSP}) {
        # + Forward strand (+ve) and primer only
        my $fsp = 0;
        my $strandc = 0;
        my $count_t = $count;
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
            # Get the suffix.
			my ($suf) = $name =~ /\.(...)/;
            if(!exists $$strand{$name}) {$strandc++;}
            if($suf && exists $$chem{primer}->{$suf}) {$fsp++;}
        }
        if($count_t == 0) {$count_t=1; }
        if(!exists $self->{fspo} && ($fsp == $count_t) && ($strandc == $count_t)) {
            if($start == $last) {
                $self->gap_format($contig,$start,"Forward strand and primer only",$eoc,$prob);	
                delete $self->{fspo};
            }else{
                $self->{fspo} = $start;
            }	
        }elsif(exists $self->{fspo} && ($fsp != $count_t || $strandc != $count_t)) {
            $self->gap_format($contig,$self->{fspo},"Forward strand and primer only",$start,$prob);
            delete $self->{fspo};
        }elsif(exists $self->{fspo} && $fsp == $count_t && $strandc == $count_t && $start == $last) {
            $self->gap_format($contig,$self->{fspo},"Forward strand and primer only",$eoc,$prob);
            delete $self->{fspo};
        }else{}
    }
    
    if($self->{RSP}) {
        # + Reverse strand (-ve) and primer only
        my $rsp = 0;
        my $strandc = 0;
        my $count_t = $count;
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
            # Get the suffix.
			my ($suf) = $name =~ /\.(...)/;
            if(exists $$strand{$name}) {$strandc++;}
            if($suf && exists $$chem{primer}->{$suf}) {$rsp++;}
        }
        if($count_t == 0) {$count_t=1; }
        if(!exists $self->{rspo} && ($rsp == $count_t) && ($strandc == $count_t)) {
            if($start == $last) {
                $self->gap_format($contig,$start,"Reverse strand and primer only",$eoc,$prob);
                delete $self->{rspo}; 
            }else{
                $self->{rspo} = $start;
            }	
        }elsif(exists $self->{rspo} && ($rsp != $count_t || $strandc != $count_t)) {
            $self->gap_format($contig,$self->{rspo},"Reverse strand and primer only",$start,$prob);
            delete $self->{rspo};
        }elsif(exists $self->{rspo} && $rsp == $count_t && $strandc == $count_t && $start == $last) {
            $self->gap_format($contig,$self->{rspo},"Reverse strand and primer only",$eoc,$prob);
            delete $self->{rspo}; 
        }else{}
    }

    if($self->{SS}) {
        # + Single Stranded regions, either forward or reverse.
        my $rev_strandc  = 0;
        my $fwd_strandc  = 0;
        my $type;
        my $count_t = $count;
	my $wgs_set = undef;
	if($self->{CMODE}) {
            foreach my $name (keys %$rlist) {
                if(exists $$wgs{$name}) {
                    $wgs_set++;
                }
            }
        } 
        foreach my $name (keys %$rlist) {
            if(exists $$wgs{$name}) {
                if($count > 1) {$count_t -= 1;}
                next;
            }
	    if($self->{CMODE} && $wgs_set) {
                #don't incrememnt.
            }else{
		if(exists $$strand{$name}) {$rev_strandc++;}
		if(!exists $$strand{$name}) {$fwd_strandc++;}
	    }
        }      
        if($count_t == 0) {$count_t=1; }
        if($fwd_strandc == $count_t) { $type = 'Forward strand only'; }
        if($rev_strandc == $count_t) { $type = 'Reverse strand only'; }      
        if(!exists $self->{singlestrand} && ($fwd_strandc == $count_t || $rev_strandc == $count_t)) {
            if($start == $last) {
                $self->gap_format($contig,$start,($fwd_strandc==$count_t)?"Forward strand only":"Reverse strand only",$eoc,$prob);
                delete $self->{singlestrand};
            }else{
                $self->{singlestrand} = $start;
                $self->{last_strand_type} = ($fwd_strandc==$count_t) ? "Forward strand only" : "Reverse strand only"; 
            }	
        }elsif(exists $self->{singlestrand} && $fwd_strandc != $count_t && $rev_strandc != $count_t) {
            $self->gap_format($contig,$self->{singlestrand},$self->{last_strand_type},$start,$prob);
            delete $self->{singlestrand};
            delete $self->{last_strand_type};
        }elsif(exists $self->{singlestrand} && ($type ne $self->{last_strand_type})){
            $self->gap_format($contig,$self->{singlestrand},$self->{last_strand_type},$start,$prob);
            $self->{singlestrand} = $start;
            $self->{last_strand_type} = $type;
        }elsif(exists $self->{singlestrand} && ($fwd_strandc == $count_t || $rev_strandc == $count_t) && $start == $last) {
            $self->gap_format($contig,$self->{singlestrand},$self->{last_strand_type},$eoc,$prob);
            delete $self->{singlestrand};
            delete $self->{last_strand_type};
        }else{}
    }
}

sub parent_over_region {
    my($rlist, $parentname) = @_;
    my $parent_count=0;
    foreach my $readname (keys %$rlist) {
        if($readname =~ /^$parentname\..*$/) {
            ++$parent_count;
        }
    }
    return $parent_count;
}

sub screen_low_quality {
  my($self,$contig,$arr,$dna,$bq) = @_;
 
  my @qual = split(" ", $$bq);                                                                                    
  for (my $i = 0; $i<scalar(@qual);$i++) {
      if($qual[$i] < $self->{Q}) {
          push @$arr,[$contig,$i+1,$i+1,"Quality value ($qual[$i]) less than $self->{Q}",$i+1];
      }
  }
}

sub screen_low_quality_wgs_data {
    my($self,$contig,$assembled) = @_;

    my (@wgs,%low_qual);
    foreach my $assembled_from (split "\n", $$assembled) {       
	next unless ($assembled_from =~ /Assembled_from/);
	my ($rd,$start,$end) = (split(" ",($assembled_from)))[1,2,3];
	next unless (($rd =~ /c\d+\.Contig.*$/) || ($rd =~ /^contig_\d+/));
	push @wgs,[$rd,$start,$end];
    }    

    foreach my $rd (@wgs) {
	my $ent = $self->{CAF}->getEntry($rd->[0]);
	my ($id) = $ent->{Sequence} =~ /Staden_id\s(\S+)/;
	my($cl,$cr) = $ent->{Sequence} =~ /Clipping QUAL (\d+) (\d+)/;
	my $dna = substr($ent->{DNA},$cl-1,$cr-$cl+1);
	my @qual = split(" ",$ent->{BaseQuality});
        # nqs pack qual to one string.
	my $q = pack("C*",map{$_+=041} @qual);
        # ...and then unpack the substring needed...
	@qual = map{unpack("C*",$_)-041} (split(//,substr($q,$cl-1,$cr-$cl+1)));	
	my $comp = ($rd->[1] > $rd->[2]);
	if($comp) {
	    $dna = rev_comp($dna);
	    @qual = reverse @qual;
	    ($rd->[1],$rd->[2]) = ($rd->[2],$rd->[1]);
	}	        
	for (my $i = 0; $i<scalar(@qual);$i++) {
	    if($qual[$i] < $self->{G}) {
		$low_qual{$rd->[1]+$i+1} = [$id,$i+1,"MouseWGS quality ($qual[$i]) less than $self->{G}"];		
	    }    			
	}
    }
    return \%low_qual;
}

sub rev_comp {
  my ($dna) = @_;
  
  $dna =~ tr/acgtACGT/tgcaTGCA/;
  return reverse($dna);
}

#
# Single base pair run search. (usually 10 or more)
#
sub single_base_pair_runs {
    my($self,$arr,$dna,$contig) = @_;
    
    $$dna =~ s/\n//g;
    foreach my $base ('A','C','G','T') {
        my $big_len = length($$dna);
        my $sub_len = $self->{B};
        my $last_i = $big_len - $sub_len;
        
        for ( my $i = 0; $i <= $last_i; $i++ ) {
            # loop as much much as we can while the base is still of the same type...
            # ..remembering to skip over pads '-'.
            my  $match = 0;
            my $j = 0;
            if($base eq substr( $$dna, $i + $j, 1 )) {
                do {
                    ++$match if ($base eq substr( $$dna, $i + $j, 1 ));                    
                    $j++;
                }while($base eq substr( $$dna, $i + $j, 1 ) 
                       || '-' eq substr( $$dna, $i + $j, 1 )); 
                
                # loose end pad if applicable.
                $j-- if ('-' eq substr( $$dna, $i + $j-1, 1 )); 
                
                if ($match >= $sub_len) {
                    my $incr = 1;
                    push @$arr,[$contig,$i + $incr,$i+$j,"Single base pair run of $match $base\'s",$i+$incr];
                    $i +=  $j;
                }
            }
        }
    }    
}

#
# TA repeat screen, more that 5 pairs, skip pads
# maybe use index. working version first..
#
sub ta_repeat_screen {
    my($self,$arr,$dna,$contig) = @_;
    
    my ( $i, $j, $match_j,$pads );
    my $big_len = length($$dna);
    my $pair = 0;
    my $valid_base = 0;
    for ( $i = 0; $i <= $big_len; $i++ ) {
        # ..remembering to skip over pads '-'.                
        for ( $j = 0, $match_j = 0; ('T' eq substr( $$dna, $i + $j, 1 ) && 'A' eq substr( $$dna, $i + $j + 1, 1 )) || 
              ('T' eq substr( $$dna, $i + $j, 1 ) && '-' eq substr( $$dna, $i + $j + 1, 1 )) || 
              $pair && '-' eq substr( $$dna, $i + $j, 1 ); $j++ ) {                        

            # skip pads if needed.
            while('-' eq substr( $$dna, $i + $j, 1 )) {
                ++$j;
            }            
            # check if we need a T
            if('T' eq substr( $$dna, $i + $j, 1 ) && !$valid_base) { 
                $valid_base++;                  
               
            # otherwise ignore and start again.
            }elsif('A' ne substr( $$dna, $i + $j, 1 ) && $valid_base) {
                $valid_base = 0;
                last;
            }
            # have a T..look ahead and skip pads.
            if(substr( $$dna, $i + $j + 1, 1 ) eq '-') {
                do{
                    ++$j;    
                }while(substr( $$dna, $i + $j + 1, 1 ) eq '-');
                # check if first base was T and next is A...
                if(substr( $$dna, $i + $j+1, 1 ) eq 'A' && $valid_base--) {
                    $pair++;      # count pairs.                             
                    $match_j = $j+1; 
                    $j++;               
                }else{                
                    $valid_base = 0;
                    last;
                }
                # ..next base is A
            }elsif(substr( $$dna, $i + $j + 1, 1 ) eq 'A' && $valid_base--) {
                $pair++;                                   
                $match_j = $j+1;  
                $j++;           
            }
        }        
        if($pair >= 5) {
            my $section = substr( $$dna, $i, $match_j );
            ($pads) = $section =~ s/-//g; # count pads...
            $j -= $pads; # ...remove from total base count.
            my $incr = 1;
            my $start = $i + $incr;
            my $end = $i+$match_j+1;
            my $bplen = $pair * 2;
            push @$arr,[$contig,$start,$end,"TA repeat (${bplen}bp)",$start];      
        }
        $i +=  $match_j;
        $pair=0;                 
    }
}


sub init_suffix_types {
    my ($self,$origin,$addi,$pcrs,$parent_sil) = @_;

  #
  # use seqchemistry table and dyeset;
  # N.B could use hardcoded hash list. 
  #
  print STDERR "\nFetching chemistry types from oracle...";
  my ($dbh,%chemlist,%plates);
  eval {
    $dbh = WrapDBI->connect('reports',{RaiseError => 1});
  };
  if($@){
    die "Unable to get handle to oracle database.$@\n";
  }	
  my $chemlist = $dbh->fetch_all(qq[
				    SELECT SUFFIX,IS_PRIMER,ID_DYESET 
				    FROM SEQCHEMISTRY
				   ]);
  
    if($self->{SIL}) {
        
        
        my $sil = $dbh->selectall_arrayref(qq[
                                              SELECT distinct A.REQUEST
                                              , A.SOURCE||(select maprow||mapcol from rowcol_map where id_map = A.MAP)
                                              SOURCENAME
                                              , NVL((SELECT bb.NAME||rcm.maprow||rcm.mapcol FROM BECKMAN_BOX bb, bb_map
                                                     bbm2, rowcol_map rcm WHERE bbm2.id_beckmanbox = bb.id_beckmanbox and
                                                     bbm2.id_map = rcm.id_map and bb.IS_ORIGINAL = 1 AND bbm2.ID_BBMAP
                                                     IN (select min(id_bbmap) from bb_map where id_subclone =
                                                         bbm.id_subclone)),A.SOURCE||(select maprow||mapcol from rowcol_map where
                                                                                      id_map = A.MAP)) BOXNAME,
                                              A.TYPE
                                              FROM bb_map bbm,(
                                                               select sr.id_request REQUEST
                                                               , NVL(sd.platename,'PCR_'||sd.clonename||'_'||sp.FWD_PRIMER||'_'||sp.REV_PRIMER)
                                                               SOURCE
                                                               , NVL(sd.id_map,0) MAP
                                                               , NVL(bb.id_beckmanbox,0) BBID
                                                               , std.description TYPE
                                                               from sil_Request sr, sil_dnasource sd, sil_pcr sp, beckman_box bb, siltypedict std
                                                               where sr.id_source = sd.id_source
                                                               and sr.id_request = sp.id_request (+)
                                                               and sd.clonename  = ?
                                                               and sd.platename = bb.name (+)
                                                               and sr.sil_type = std.id_dict
                                                               )A
                                              where A.BBID = bbm.id_beckmanbox (+)
                                              and bbm.id_map (+)= A.MAP                                              
                                              ],{},$self->{PROJECT});


        my $sil2 = $dbh->selectall_arrayref(qq[
                                                  select rr.dnaplate_name||rm.maprow||rm.mapcol ORIGINAL,
                                                  bb.name||rm2.maprow||rm2.mapcol RETRANS
                                                  from retrans_request rr,
                                                  rrplate_bb rb,
                                                  beckman_box bb,
                                                  bb_map bm,
                                                  subclone s,
                                                  ligation l,
                                                  rowcol_map rm,
                                                  rowcol_map rm2
                                                  where rr.dp_idmap = rm.id_map
                                                  and rr.rp_idmap = rm2.id_map
                                                  and rr.id_rrplate = rb.id_rrplate
                                                  and rb.id_beckmanbox = bb.id_beckmanbox
                                                  and bb.is_original = 0
                                                  and bb.id_beckmanbox = bm.id_beckmanbox
                                                  and rr.rp_idmap = bm.id_map
                                                  and bm.id_subclone = s.id_subclone
                                                  and s.id_ligation = l.id_ligation
                                                  and l.clonename = ?
                                                  order by bb.name, rm2.maprow, rm2.mapcol
                                                  ],{},$self->{PROJECT});


        my $retrans =  $dbh->prepare(qq[
                                        select rr.dnaplate_name||rm.maprow||rm.mapcol ORIGINAL,
                                        bb.name||rm2.maprow||rm2.mapcol RETRANS
                                        from retrans_request rr,
                                        rrplate_bb rb,
                                        beckman_box bb,
                                        bb_map bm,
                                        subclone s,
                                        ligation l,
                                        rowcol_map rm,
                                        rowcol_map rm2
                                        where rr.dp_idmap = rm.id_map
                                        and rr.rp_idmap = rm2.id_map
                                        and rr.id_rrplate = rb.id_rrplate
                                        and rb.id_beckmanbox = bb.id_beckmanbox
                                        and bb.is_original = 0
                                        and bb.id_beckmanbox = bm.id_beckmanbox
                                        and rr.rp_idmap = bm.id_map
                                        and bm.id_subclone = s.id_subclone
                                        and s.id_ligation = l.id_ligation
                                        and l.clonename = ?
                                        and bb.name = ?
                                        order by bb.name, rm2.maprow, rm2.mapcol                                        
                                        ]);
        # Pain and Suffering Pt.2.
        # Stores SILs in hash but assigns each read to its parent.
        # fix_read everywhere as the DB may not hold what is actually on disk. (!expletive!)
        foreach my $ref (@$sil2) {
            my($orig,$retrans) = @$ref;
            $self->fix_read(\$orig);
            $self->fix_read(\$retrans);
            $$origin{$retrans} = $orig;
        }

        my %tocheck;
        foreach my $ref (@$sil) {
            my($request,$source,$parent, $type) = @$ref;
            # skip restriction fragment sils. No need to tag.
            next if ($type eq 'Restriction');
            # normal sil with parent.
            if($source ne $parent) {
                $self->fix_read(\$request);
                $self->fix_read(\$parent);
                $self->fix_read(\$source);
                $$origin{$request} = $parent;
                $$origin{$source}  = $parent;
            }else{                
                # re-iterate.
                $self->fix_read(\$parent);
                $tocheck{$request} = $parent; # store for checking.
                # job on another job.
                if($parent =~ /^(R\d+)[a-p]\d{1,2}.*$/) {
                    $retrans->execute($self->{PROJECT},$1);
                    my ($info) = $retrans->fetchall_arrayref();                  
                    if(defined $info) {
                        foreach my $rt (@{$info}) {
                            $$origin{$rt->[1]} = $rt->[0];                            
                            $self->fix_read(\$rt->[0]);
                            $self->fix_read(\$rt->[1]);          
                        }                                            
                    }                   
                }else{
                    # must be pcr for sil on clone dna.
                    $self->fix_read(\$request);
                    $self->fix_read(\$parent);
                    $$origin{$request} = $parent;
                    # add for pcr as is pcr_for_sil
                    $$pcrs{$request}++;
                    
                }
            }
        }
        $retrans->finish(); 

        # HACK. for checking those jobs that dont fall into queries above
        # but can work out info from whats already found. see zK6H12 or kva
        
        # occurs where one entry is read  - parent and another is
        #                           read2 - read
        # so we iterate over and assign read2 - parent.
        while(my($job,$parent) = each %tocheck) {
            while(my($key,$val)= each %$origin) {
                if($val =~ /^\D(\d+(\w))?(\-)?(\S+)?[a-p]\d{1,2}..*$/) { 
                    if($1 == $job) {
                        $$origin{$job} = $$origin{$parent};
                        $$origin{$key} = $$origin{$parent};
                    }
                }
            }
        }



    }

  
  if($self->{ADDI} ||$self->{UADDI} )  {
      
      my $additive = $dbh->prepare(qq[
                                      select a.description
                                      from seqread sr, gel_sample gs, request_map rm, finish_request fr,
                                      additive a
                                      where sr.name = ?
                                      and a.id_additive in (1,2,3,4,14,15,16,25,29,30,31)
                                      and sr.id_gel = gs.id_gel
                                      and sr.lane = gs.lane
                                      and gs.id_loadingsample = rm.id_loadingsample
                                      and rm.id_frequest = fr.id_frequest
                                      and fr.id_additive = a.id_additive
                                      ]);
      
      foreach my $plate ($self->{CAF}->readsList()) {
	  $additive->execute($plate);
	  my ($addi_desc) = $additive->fetchrow_array();
	  if(defined $addi_desc) {
	      $$addi{$plate}++;  
          }
      }
      $additive->finish();
  }        
  
   if($self->{PCR}) {
       my $pcr = $dbh->selectall_arrayref(qq[
                                             select pcr_template_name from
                                             pcr_product where dna_source = ?
                                             ],{},$self->{PROJECT});

       
       map { $$pcrs{$_->[0]}++ } @$pcr;
   }

  print STDERR "Done\nBuilding lookups...";

  # Build a hash containing all reads that are 
  # terminator and of those, all reads that are 
  # big dyes. Anything not in the list is taken
  # as a primer.
  foreach my $cline (@$chemlist) {
      if($cline->[1] == 0) {
          $chemlist{term}->{$cline->[0]}++;
          # 3 for normal BDT, 7 for BDT VIII.
          if(($cline->[2] == 3) || ($cline->[2] == 7)) { 
              $chemlist{bdt}->{$cline->[0]}++;
          }    
      }else{
          $chemlist{primer}->{$cline->[0]}++;
      }
  }
  
  print STDERR "Done\n";
  $dbh->disconnect();

  return \%chemlist;
  
}

sub parent_sil {
    my($self,$silhash,$name,$job) = @_;
    
    my $orig = $name;
    if($orig =~ /(\S+[a-p])0(\d{1}..*)$/) { 
        $orig = $1.$2;
    }    
    if($name =~ /^$$silhash{$job}\..*$/) {	    
        #exists as a parent read.
        return 1;
    }elsif($orig =~ /^$$silhash{$job}\..*$/) {
        #..try without the leading zero...
        return 1;
    }else{
        return 0;
    }
}

sub fetch_region {
  my ($self,$pos,$contig_reads,$reg) = @_;

  #
  # Find region, looping to allow multiple reads starting
  # and ending at the same point.

  # @$contig_reads is the array generated from the get_assembled_from routine below.

  do{
    my $data = shift(@$contig_reads);  
    # Delete the read if it exists in the region hash.
    if (exists $$reg{$data->[0]}) {
      delete $$reg{$data->[0]};
      $$pos = $data->[1];
    }else{    
      # If it's a new read, add to the hash and 
      # alter the new position.
      $$reg{$data->[0]}++;      
      $$pos = $data->[1];
    }    
  }while(@$contig_reads && $$contig_reads[0]->[1] == $$pos);
  return;
}

sub get_assembled_from {
    my ($self,$seq) = @_;
    
    my $av = -1;
    my (@af,$last,%strand);
    
    foreach my $line (split("\n", $$seq)) {    
        my ($af, $read, $cs, $ce, $rs, $re) = split(/\s+/, $line);
        next unless ($af eq "Assembled_from");    
        my $complement = ($cs > $ce);    
        if ($complement) {
            ($cs, $ce) = ($ce, $cs);
            $strand{$read}++ if $complement
            }

        # Add the start of a read and the end seperatly for sorting by position later.        
        push(@af,[$read,$cs-1]);  
        push(@af,[$read,$ce]);
    }  
    
    unless (@af) { return ([], {}, 0, 0); }  # Nothing in this contig.
    
    @af = sort along_contig @af;  
    
    # Find last postion start and end of contig.
    # Loop to allow for reads ending at the same point.
    
    for(my $av=@af-1;$av>=0;$av--) {
        $last = $af[$av]->[1];
        last if ($last != $af[$#af]->[1]);
    }

    return (\@af,\%strand,$last,$af[$#af]->[1]);
}

sub along_contig {
  my ($self) = shift;
  return ($a->[1] <=> $b->[1]);
  
}


sub parse_args {
  my($self,$set) = @_;

  if ($set eq 'all') {
      # Make all searches available.
      foreach my $search (qw(ST 2R 3R ADDI UADDI SIL PCR RR+1 WW+1 WGS+1 FSP RSP FST RST QUAL+30 TAREP SBPR+10)) {
          if ($search =~ /RR/) { ($self->{RR},$self->{R}) = split/\+/,$search;}
          elsif ($search =~ /WW/) { ($self->{WW},$self->{W}) = split/\+/,$search;}
          elsif ($search =~ /WGS/) { ($self->{WGS},$self->{M}) = split/\+/,$search;}
          elsif ($search =~ /QUAL/) { ($self->{QUAL},$self->{Q}) = split/\+/,$search;}
          elsif ($search =~ /SBPR/) { ($self->{SBPR},$self->{B}) = split/\+/,$search;}	  
          else { $self->{$search}++; }
      }
  }else {
      # Use the searches selected.
      foreach my $search (split/,/,$set) {        
          if ($search =~ /RR/) {($self->{RR},$self->{R}) = split/\+/,$search;}
          elsif ($search =~ /WW/) {($self->{WW},$self->{W}) = split/\+/,$search;}
          elsif ($search =~ /WGS\+/) {($self->{WGS},$self->{M}) = split/\+/,$search;}
          elsif ($search =~ /QUAL/) { ($self->{QUAL},$self->{Q}) = split/\+/,$search;}
          elsif ($search =~ /SBPR/) { ($self->{SBPR},$self->{B}) = split/\+/,$search;}
	  elsif ($search =~ /WGSS\+/) { ($self->{WGSS},$self->{G}) = split/\+/,$search;}
          else { $self->{$search}++; }      
      }
  }
}

sub gap_format {
  my($self,$ct,$start,$string,$to,$ray) = @_;

  if($string eq 'SIL:Single') {
      push @$ray, [$ct,$start+1,$to,'SIL only',$start+1];
      push @$ray, [$ct,$start+1,$to,'Single template',$start+1];
  }else{
      push @$ray, [$ct,$start+1,$to,$string,$start+1];
  }
}


sub sort_list {
  my ($self,$list) = @_;
  
  my @sorts;
 
  # parse and join consecutive problem types. used for SIL:Single.
   #see above.
  my %types;
  map { push @{$types{$_->[3]} }, $_ } @$list;
  
  foreach my $prob (keys %types) {
      my $last_start;
      my $last_ctg;
      my $index=0;
#     foreach my $prob_ref ( sort{ $a->[1] <=> $b->[1] } @{$types{$prob}} ){          
      my @js10a = sort{ $a->[1] <=> $b->[1] } @{$types{$prob}};
      foreach my $prob_ref ( @js10a ){          
          if(!defined $last_start || $prob_ref->[1] != $last_start+1 && $prob_ref->[0] eq $last_ctg) { 
              $last_start = $prob_ref->[2];              
              $last_ctg = $prob_ref->[0];
              ++$index;
              next; 
          }else{ # join segment.
              my $first = @{$types{$prob}}->[$index-1];
              my $curr_end = $prob_ref->[2];
              $first->[2] = $curr_end;
              splice @{$types{$prob}}, $index-1, 2, $first;              
              $last_start = $curr_end;
              $last_ctg = $prob_ref->[0];
              # no need to increment index as we've spliced the array.
          }
      }
  }
  
  my @joined;
  map { map { push @joined, $_ } @{$types{$_}} } keys %types;
    
  # quick map-sort-map.
  push @sorts, map{ $_->[0] } sort{ $a->[1] <=> $b->[1] } map{[ $_, $_->[4] ]} @joined;

  # Assign our sorted list to the original.
  @$list = @sorts;

}



1;

=head1 NAME

 GNUtils - GapNav utility module.

=head1 SYNOPSIS

 use GNUtils;
 use Caftools;

 my $caf = Caftools->new();
 $caf->readCafFile($file);

 my $gnav = GNUtils->new($caf);
 my $problems_ref = $gnav->find_problem_regions(OPTS);

 e.g: 
  my $problems_ref = $gnav->find_problem_regions(ST,RR+1,QUAL+40)
  
  my $problems_ref = $gnav->find_problem_regions($caf, all);


  $gnav->Visualize($problems,"mygif.gif");


=head1 OPTIONS

=over 4

=item all    = Overides all other options and sets the args to default.
 (ST 2R 3R SS ADDI SIL PCR RR+1 WW+1 FST RST FSP RSP QUAL+30)

=item ST     = Single template regions

=item 2R     = Two read only regions

=item 3R     = Three read only regions

=item SS     = Single stranded regions

=item ADDI   = Additive only regions (tagged by additive2gapdb)

=item UADDI  = Uni-directional Additive only regions (tagged by additive2gapdb or asp)

=item SIL    = SIL only regions with below selections. Stores type (TIL/SIL) also works
               on whether the reads have parent subclones over region.

=item PCR    = PCR only regions

=item RR     = ReadRaid only or (RR+1) ReadRaid + 1 project read.

=item WW     = WWRaid only or (WW+1) WWRaid + 1 project read.

=item FST    = Forward strand and terminator only

=item RST    = Reverse strand and terminator only

=item FSP    = Forward strand and primer only

=item RSP    = Reverse strand and primer only

=item QUAL   = Search for consensus bases < Q 

=item SBPR   = Search for single base pair runs of and base.

=item TAREP  = Search for TA repeats. usaully more that 5 pairs long.

=item WGSS   = WGS quality search.

=back

=head1 DESCRIPTION

 Utility module mainly used by the program 'gapNav'. It will take a 
 caf object and run a number of user selected searches on that caf file
 by splitting the contigs up into read regions. This is best defined as 
 an area between 1 reads end and another reads start. 

 It returns a reference to an array or arrays containing the positions found in
 the format...

  [ contig, start, end, type (descriptive string) ]

  e.g   [00175,1,706,ReadRaid only] 

 This array can be formatted and written to a file for navigation 'by file' within 
 the search menu of gap4.

 There is also the Visualize routine which is passed the problems list
 and the name of an image to write. This will generate a representation 
 of the database with problem regions plotted.


 Each type of problem search follows the same pattern. e.g Single Template

    if($self->{ST}) {
        # Single Template.
        my %seen=();
        foreach my $name (keys %$rlist) {
            # cope with leading zeros etc...
            $self->fix_read(\$name);
            my ($tp,$sf) = split/\./,$name;
            $seen{$tp}++;      
        }       
        $self->{current_t} = keys %seen;
        # in problem region, start not set.....
        if(!exists $self->{stset} && $self->{current_t} == 1) {
            # print if end of contig.
            if($start == $last){
                $self->gap_format($contig,$start,"Single template",$eoc,$prob); 
                delete $self->{stset};
            }else{
                # problem region....set start....
                $self->{stset} = $start;
            }              
            # start set but out of problem region.... so report..
        }elsif(exists $self->{stset} && $self->{current_t} != 1) {
            $self->gap_format($contig,$self->{stset},"Single template",$start,$prob);
            delete $self->{stset};
            
            # start set, still in problem region but at end of contig.... so report...
        }elsif(exists $self->{stset} && $self->{current_t} == 1 && $start == $last) {
            $self->gap_format($contig,$self->{stset},"Single template",$eoc,$prob);
            delete $self->{stset};
            # do nothing ... keep iterating ....
        }else{}
    }



=head1 AUTHOR

 Mark Griffiths

=head1 HISTORY

 #
 # $Log: GNUtils.pm,v $
 # Revision 1.114  2006/10/25  15:54:19  js10
 # Marks changes
 #
 # Revision 1.114  2006/10/24 14:10:59  js10
 # Finished changes started by mng
 #
 # Revision 1.112  2005/11/18  14:54:19  mng
 # Added TA repeat search
 #
 # Revision 1.73  2004/02/19  11:19:51  mng
 # bug fix on error warnings for sils.
 #
 # Revision 1.72  2003/12/10  08:19:30  mng
 # Added sql stuff for finding new style SILS from oracle.
 #
 # Revision 1.71  2003/05/14  08:13:54  mng
 # Fixed problem where GD has removed gif support in favour for png.
 #
 # Revision 1.70  2003/05/13  08:18:57  bt1
 # *** empty log message ***
 #
 # Revision 1.69  2003/02/06 13:01:43  mng
 # Fixed bug in clipping segment for WGS quality search.
 #
 # Revision 1.68  2003/02/06  08:16:48  mng
 # *** empty log message ***
 #
 # Revision 1.67  2003/02/04  14:14:42  mng
 # *** empty log message ***
 #
 # Revision 1.66  2003/02/03  11:35:22  mng
 # Pattern match fix for wgs low qual regions.
 #
 # Revision 1.65  2003/01/30  13:38:20  mng
 # *** empty log message ***
 #
 # Revision 1.64  2003/01/29  13:33:06  mng
 # *** empty log message ***
 #
 # Revision 1.63  2003/01/29  13:31:10  mng
 # *** empty log message ***
 #
 # Revision 1.62  2003/01/29  13:21:18  mng
 # Added WGS search for quality less than user-specified value.
 #
 # Revision 1.61  2003/01/15  09:53:48  mng
 # removed term finding throughout the database if in checking mode.
 #
 # Revision 1.60  2003/01/14  10:48:43  mng
 # Updated to remove single stranded region listings from those areas
 # covered by wgs data.
 #
 # Revision 1.59  2002/10/07  13:19:11  mng
 # Fixed a couple of warning messages.
 #
 # Revision 1.58  2002/10/07  06:57:12  mng
 # *** empty log message ***
 #
 # Revision 1.57  2002/10/07  06:42:59  mng
 # *** empty log message ***
 #
 # Revision 1.56  2002/10/04  08:18:35  mng
 # Added code for checking mode. skips regions covered by wgs and fsbdt/rsbdt.
 #
 # Revision 1.55  2002/09/19  09:53:51  mng
 # Fixed bug in search routines where a wgs only region would show up as
 # all problems in the tag section. Bug due to wgs ignored so routines
 # would succeed as 0 reads = 0 problems.
 #
 # Revision 1.54  2002/08/22  14:22:00  mng
 # Reworked single base pair run routine so that it can cope
 # with any number of pads.
 #
 # Revision 1.53  2002/08/06  14:16:25  mng
 # Added/altered code to allow selection of readraid/wwraid and wgs regions with ?
 # templates defined on main interface.
 #
 # Revision 1.52  2002/07/10  06:58:20  mng
 # Fixed pattern match for SIL job
 #
 # Revision 1.51  2002/07/09  09:11:05  mng
 # Fixed but for WGS data. Always set read count to 1 so that wgs
 # single template regions don't become zero and match everything.
 #
 # Revision 1.50  2002/07/08  07:38:12  mng
 # Updated to remove wgs data.
 #
 # Revision 1.43  2002/06/11  06:59:54  mng
 # *** empty log message ***
 #
 # Revision 1.42  2002/06/11  06:33:51  mng
 # padded bug fix on snbp routine.
 #
 # Revision 1.41  2002/06/10  14:39:50  mng
 # Re-worked string searching algorithm to take pads into account.
 #
 # Revision 1.40  2002/06/10  08:59:22  mng
 # Added single base pair run checks. Doesn't use -cutoff Q to translate
 # and index on n's for quality either.
 #
 # Revision 1.39  2002/04/22  11:41:44  mng
 # *** empty log message ***
 #
 # Revision 1.38  2002/04/18  14:49:18  mng
 # Added fetchwgs readhash entries for type FWGS
 #
 # Revision 1.37  2002/04/17  13:34:18  mng
 # Fixed WGS bug
 #
 # Revision 1.36  2002/04/05  07:14:23  mng
 # Added arachne assembly naming convention (contig_1234).
 #
 # Revision 1.35  2002/03/04  09:22:26  mng
 # Added functions to find MouseWGS data + project reads.
 #
 # Revision 1.34  2002/01/28  09:04:17  mng
 # *** empty log message ***
 #
 # Revision 1.32  2002/01/09  08:07:44  mng
 # Added yet more parsing of reads to find exact SIL job. No
 # looks for J/T + numbers + letters (+possibly -0).
 #
 # Revision 1.31  2002/01/08  09:53:34  mng
 # Added BDT VIII using id_dyeset 7 within init_suffix_types.
 #
 # Revision 1.30  2001/12/12  12:20:25  mng
 # Fixed no type navigation file line for UADDI.
 #
 # Revision 1.29  2001/12/12  11:06:00  mng
 # Re-worked SIL sql. Makes finding the parent read more accurate.
 #
 # Revision 1.28  2001/12/11  11:44:39  mng
 # Added/fixed uni-directional additive + primer section.
 #
 # Revision 1.27  2001/12/06  10:13:44  mng
 # Fixed SIL parent reads and added R reads o SIL regions.
 #
 # Revision 1.26  2001/11/16  11:08:23  mng
 # *** empty log message ***
 #
 # Revision 1.25  2001/11/16  07:49:34  mng
 # Added UADDI code to find uni-directional additives.
 #
 # Revision 1.24  2001/10/25  15:48:40  mng
 # single strand option is only available if you specifically
 # specify it otherwise it may cause problems.
 #
 # Revision 1.23  2001/10/25  12:25:47  mng
 # Added single strand option.
 #
 # Revision 1.22  2001/10/04  10:59:02  mng
 # Fixed consensus ignoring bug.
 #
 # Revision 1.21  2001/10/02  06:34:38  mng
 # Altered ADDI code to find additive reads from oracle...this is a bit slow!!
 # Now include additive reads in primer only regions.
 #
 # Revision 1.20  2001/09/26  07:31:00  mng
 # Fixed SIL problems with wells
 #
 # Revision 1.19  2001/09/26  07:14:02  mng
 # Added real read to sil_parent has as well as that with a
 # 0 prepended to wells 1-9
 #
 # Revision 1.18  2001/09/18  13:26:43  mng
 # Now ignores consensus reads.
 #
 # Revision 1.17  2001/09/13  10:24:57  mng
 # Added T reads for sil_parent searches.
 #
 # Revision 1.16  2001/09/05  14:09:30  mng
 # Fixed primer hash tables.
 #
 # Revision 1.15  2001/09/04  10:55:11  mng
 # Fixed parent_sil pattern match as it needed to be anchored at the front
 # 'cuz it was matching any template with the correct pattern.
 #
 # Revision 1.14  2001/08/15  07:05:27  mng
 # Added function to find SIL parent and include that in a SIL only region
 # 
 #
 # Revision 1.13  2001/07/31  14:22:00  mng
 # Fixed code parsing bug
 #
 # Revision 1.12  2001/07/30  14:20:39  mng
 # Fixed sorting bug and init variables errors.
 #
 # Revision 1.11  2001/07/10  14:31:31  mng
 # *** empty log message ***
 #
 # Revision 1.10  2001/07/03  06:32:59  mng
 # Fixed usless variable warnings.
 #
 # Revision 1.9  2001/06/28  14:36:45  rmd
 # Removed debugging print statements.
 #
 # Revision 1.8  2001/06/28  14:30:24  rmd
 # Fixed problem where get_assembled_from died if there were no assembled_from
 # lines.
 #
 # Revision 1.7  2001/06/28  13:55:13  mng
 # Changed while loop to find last and end positions in contig
 # for 'for' loop, added revision logs to pod history.
 #
 # Revision 1.6  2001/06/28  13:45:29  mng
 # *** empty log message ***
 #
 # Revision 1.5  2001/06/28  11:23:10  mng
 # bug fixes
 #
 # Revision 1.4  2001/06/27  14:29:58  mng
 # *** empty log message ***
 #
 # Revision 1.3  2001/06/27  14:27:51  mng
 # *** empty log message ***
 #
 #  Revision 1.2  2001/06/27  14:20:48  mng
 #  changed array of strings returned to array of arrays.
 # 
 # Revision 1.1  2001/06/26  09:10:49  mng
 # Initial revision
 #
