#!/usr/bin/perl

use strict;
use warnings;

use lib '../';

#use ArbitraryCoord;
use CPDNucReads;

my $numbins = 60;
my $binwidth = 1000;
# ask for probe filename to analyze
print STDERR "Enter filename of plus strand reads\n";
my $plusfile = <STDIN>;
chomp($plusfile);

print STDERR "Enter filename of minus strand reads\n";
my $minusfile = <STDIN>;
chomp($minusfile);

# add filename here for probe values 
#my $filename = "h2bwtdata.txt";

print STDERR "Loading Probe Values\n";
my $reads = CPDNucReads->new($plusfile, $minusfile);

#print header
print "Telomere data\nSequencing data from file: $plusfile\t$minusfile\n";
print "Distance from telomere end\n";
print "Measurement";
for (my $i = 0; $i < $numbins; $i++ )
{
	my $dist = ($i * $binwidth ) + ($binwidth / 2.0);
	print "\t$dist";
	
}
print "\n";

my @cpdval;
my @cpdcount;

open ( CHROMS, "../ChromInfo.txt") || die "Couldn't open chromosome file\n";
my %chrs;
while ( <CHROMS> )
{
	chomp $_;
	my @temp = split /\t/, $_;
	if ( $temp[0] =~ /chr[XIV]+/ )
	{
		$chrs{$temp[0]} = $temp[1];
	}
}


foreach my $chr (sort keys %chrs)
{
	print STDERR "Starting $chr\n";
	my %plusreads = $reads->get_plus_reads_for_chromosome($chr);
	my $num_plusreads = scalar keys %plusreads;
	my %minusreads = $reads->get_minus_reads_for_chromosome($chr);
	my $num_minusreads = scalar keys %minusreads;
	print STDERR "$chr reads: $num_plusreads plus reads and $num_minusreads minus reads\n";

	for (my $i = 0; $i < $numbins; $i++ )
	{
		# calculate left telomere
		my $start = $i * $binwidth + 1;
		my $end = ($i + 1) * $binwidth;
		for ( my $pos = $start; $pos <= $end; $pos++ )
		{
			if ( exists $plusreads{$pos} )
			{
				$cpdval[$i] += $plusreads{$pos};
				$cpdcount[$i]++;
			}
			if ( exists $minusreads{$pos} )
			{
				$cpdval[$i] += $minusreads{$pos};
				$cpdcount[$i]++;
			}
		
		}
		my $chrend = $chrs{$chr};
		$start = $chrend - (($i + 1) * $binwidth) + 1;
		$end = $chrend - ($i * $binwidth);
                for ( my $pos = $start; $pos <= $end; $pos++ )
                {       
                        if ( exists $plusreads{$pos} )
                        {       
                                $cpdval[$i] += $plusreads{$pos};
                                $cpdcount[$i]++;
                        }       
                        if ( exists $minusreads{$pos} )
                        {       
                                $cpdval[$i] += $minusreads{$pos};
                                $cpdcount[$i]++;
                        }       

                }  
	}

}
print "Average:";
for (my $i = 0; $i < $numbins; $i++ )
{
	my $mean = 1.0 * $cpdval[$i] / $cpdcount[$i];
	print "\t$mean";
}
print "\nCPD count:";
for (my $i = 0; $i < $numbins; $i++ )
{
        print "\t$cpdval[$i]";
}
print "\nDIPY count:";
for (my $i = 0; $i < $numbins; $i++ )
{
        print "\t$cpdcount[$i]";
}
print "\n";


