#!/usr/bin/perl

use strict;
use Getopt::Long;
my $minlen=100;
my $thresh=1.0;
my $blatfile;

GetOptions( "blat=s"=>\$blatfile, "thresh=s"=>\$thresh, "minlen=s"=>\$minlen );

open(PSL,$blatfile) || die "Could not open $blatfile\n";

my $total = 0;
my $cover;
my $lines=0;
my $q=0;
while(<PSL>) {
  chomp;
  my @a = split;
  my ( $match, $qlen, $qfrom, $qto, $contig, $chr, $len, $from, $to ) = ($a[0], $a[10], $a[11], $a[12], $a[9], $a[13], $a[12], $a[15], $a[16] ) ;
#  print join(  "\t", $match, $contig, $chr, $len, $from, $to ) . "\n";
  my $tlen = abs($to-$from);
  if ( $match >= $thresh*$qlen && $qlen>=$minlen ) {
    $total += $len;
    if ( $from < $to ) {
      push @{$cover->{$chr}}, [$from, $to];
    }	
    else {
      push @{$cover->{$chr}}, [$to, $from];
    }	
  }
  elsif ( $qlen>500 ) {
    print "$_\n";
#    print join( "\t", $q++, $match, $qlen, $qfrom, $qto, $tlen, $contig, $chr, $len, $from, $to ) . "\n";
  }
	
  $lines++;
  print "$lines\n" if ( !($lines % 100000));
}


foreach my $chr ( sort keys %$cover ) {
  print "scanning $chr\n";
  my @c;
  my @sorted = sort { my $x = $a->[0] <=> $b->[0]; $x = $a->[1] <=> $b->[1] if ( !$x); return($x) } @{$cover->{$chr}};

  foreach my $y ( @sorted ) {
#   print "$y->[0] $y->[1]\n";
    for( my $i=$y->[0];$i<=$y->[1];$i++) {
      $c[$i]++;
    }
  }

  my $k=0;
  my $m=0;
  my $u;
  while($c[$k]==0 && $k < $#c ) {$k++};

  for( my $i=$k;$i<=$#c;$i++) {
    $m += $c[$i]>0;
    $u += $c[$i]==1;
    
  }	
  my $len = $#c-$k+1;
  my $x = $m/$len;
  my $w = $u/$len;
  print "$chr $k $#c $len $m $x $u $w\n";
}
