

use strict;
use warnings;

use Getopt::Std;

our($opt_m, $opt_l);

getopts('m:l:');

# Bd1     27845132        A       13      .............   ~~~~~~~~~~~~~



my $MIN_COVERAGE = $opt_m || 30;
my $MIN_LENGTH = $opt_l || 30;

print STDERR "Using a min coverage value of $MIN_COVERAGE\n";
my $previous_count=0;
my $previous_coord=0;
my $previous_id="";
my $start_coord = 0;
my $end_coord =0;
my $region_length = 0;
my $pileupsum =0;

while (<>) { 
    chomp();
    my ($id, $coord, $nuc, $count, $align, $qual) = split /\t/;

    if ($count >= $MIN_COVERAGE) { 
	# in a region
	$region_length++;
	$pileupsum += $count;
	if ($previous_count < $MIN_COVERAGE) { 
	    # 
	    #print STDERR "Starting new region\n";
	    $region_length =1;
	    $start_coord = $coord;
	}
	#print STDOUT "We are in a region ($id, $coord, $region_length, $pileupsum, $pileupsum, $previous_count, $start_coord)\n";
	
    }
    else { 
	
	# ending a region
	if ($start_coord) { 
	    $end_coord = $coord -1;
	    my $density = $pileupsum / $region_length;

	    if ($MIN_LENGTH < $region_length) { 
		print STDOUT "$id, $region_length, $start_coord, $end_coord, $density\n";
	    }
	    $region_length = 0;
	    $density = 0;
	    $start_coord = 0;
	    $end_coord = 0;
	    $pileupsum =0;
	
	}
	   
    }

    $previous_coord = $coord;
    $previous_count = $count;
    $previous_id = $id;

}
    
