# use strict;
#
my $f = 0;
my $ch = 0;
my $element = $ARGV[0];
my $sample = $ARGV[1];
my $inputstep = 50;
my $inputwin = 50;
my $exten_fold = 1;
my $win_star =0 ;
my $seq=0;
my @chrom = qw(chr1 chr2 chr3 chr4 chr5 chr6 chr7 chrX chr8 chr9 chr10 chr11 chr12 chr13 chr14 chr15 chr16 chr17 chr18 chrY chr19 chr20 chr21 chr22);


my $out_file = $sample. "." . "$element" . ".ud10k";
open my $outfile1, ">", "./$out_file"  or die $!; #input bowtie out file

for ($ch = 0; $ch <= $#chrom; $ch++){
$CHR = $chrom[$ch];

open my $chrfile, "<", "../reads_pos/$sample/$CHR"  or die $!; #input bowtie out file
open my $locfile, "<", "./elements/$element/$CHR" or die $!; # end,strand,pinum,clone,contigs_num(from class_10k_50pi)

my @reads_pos = ();chomp ( @reads_pos = <$chrfile> );@reads_pos = sort { $a <=> $b } @reads_pos;

while (my $line = <$locfile>){
	chomp $line;
my @arr = split /\t/, $line;
	$class_name = $arr[0];
	$class_type = $arr[3];
	$class_star = $arr[1];
  	$class_end = $arr[2];
  $ref_len = $class_end - $class_star + 1;
print $outfile1 "$class_name\t$ref_len";
#print "$class_name\t$class_type\t$ref_len\n";
my $class_len = 10000; my $all_star = $class_star - ($class_len*$exten_fold); my $all_end =  $class_end + ($class_len*$exten_fold);

my @sub_reads_pos = ();foreach $reads_pos (@reads_pos){if($reads_pos >= $all_star - 1000 and $reads_pos <= $all_end + 1000){push @sub_reads_pos, $reads_pos;}}
@sub_reads_pos = sort { $a <=> $b } @sub_reads_pos;

undef @numbers;
## left ## 50 window, 
my $seq = 0; my $density = 0; 
$left_star = $class_star - 10000;$win_star = $left_star;$w=$class_len/$inputwin;$step=$class_len/$inputstep;
for($j=0;$j<$inputstep;$j++){
$win_end = $win_star + $w;
foreach $sub_reads_pos (@sub_reads_pos){
if($sub_reads_pos >= $win_star and $sub_reads_pos < $win_end){
$seq++;
}
}
$seq_norm=$seq/$w*1000;
$seq=sprintf("%.3f","$seq_norm");
$numbers[$j] = $seq;
$seq = 0;
$win_star = $win_star + $step;
}
## gene_body
my $seq = 0; my $density = 0;
$win_star = $class_star; $w=$ref_len/$inputwin; $step=$ref_len/$inputstep;
for($j=$inputstep;$j<$inputstep*2;$j++){
$win_end = $win_star + $w;
foreach $sub_reads_pos (@sub_reads_pos){
if($sub_reads_pos >= $win_star and $sub_reads_pos < $win_end){
$seq++;
}
}
$seq_norm=$seq/$w*1000;
$seq=sprintf("%.3f","$seq_norm");
$numbers[$j] = $seq;
$seq = 0;
$win_star = $win_star + $step;
}
##right
my $seq = 0;
$win_star = $class_end;$w=$class_len/$inputwin; $step=$class_len/$inputstep;
for($j=2*$inputstep;$j<3*$inputstep;$j++){
$win_end = $win_star + $w;
foreach $sub_reads_pos (@sub_reads_pos){
if($sub_reads_pos >= $win_star and $sub_reads_pos < $win_end){
$seq++;
}
}
$seq_norm=$seq/$w*1000;
$seq=sprintf("%.3f","$seq_norm");
$numbers[$j] = $seq;
$seq = 0;
$win_star = $win_star + $step;
}

if($class_type eq "-"){@numbers = reverse @numbers;}
for($j=0;$j<3*$inputstep;$j++){
#print"\t$numbers[$j]";                 
print $outfile1 "\t", $numbers[$j];
}
print $outfile1 "\n";
undef @numbers;
#print "\n";
}
}
#close $density_out_file;
$avg_len=$inputstep*3;
system "mv $out_file ./density_files/";
#}