#!/usr/bin/perl -w 

print "prob best to do it from bed files\n";
print "doing it from bed files\n";
print "ALL POS STRAND\n";
#has to be ints
$window = 30;
$half = 15;
$twice = 60;
print "WINDOW $window\n";

@files = qw(sr_2AB_1_BC1_sorted.bed
sr_2AB_2_BC2_sorted.bed
sr_2AB_T_BC3_sorted.bed
sr_45_2_BC5_sorted.bed
sr_45a_1_BC4_sorted.bed
sr_45b_1_BC6_sorted.bed
sr_91a_1_BC7_sorted.bed
sr_91a_2_BC8_sorted.bed
sr_91b_1_BC9_sorted.bed
sr_91b_2_BC10_sorted.bed);


$color[0]= "150,50,50";
$color[1]= "50,150,50";
$color[2]= "50,50,150";
$color[3]= "150,150,50";
$color[4]= "150,50,150";
$color[5]= "50,150,150";
$color[6]= "200,30,30";
$color[7]= "30,200,30";
$color[8]= "30,30,200";
$color[9]= "50,250,250";


open(TK,">allfractracks.txt");
#print TK "track all_spljun_frac
#type bigWig 
#container multiWig
#aggregate transparentOverlay
#showSubtrackColorOnUi on
#viewLimits 0:1
#";

$nn=0;
foreach $f (sort @files){
%total = ();
%spljun = ();

#parent all_spljun_frac
$short = $f;
$short =~ s/_sorted.bed//;
$short = "_frac.".$short;

print TK "
track $short
shortLabel $short
longLabel $short
type bigWig 
bigDataUrl http://intron.ucsc.edu/jen/may2021/in_frac.$f.bw
visibility full
alwaysZero on
viewLimits 0:1
color $color[$nn]

";
$nn++;

$max=0;$min=40000;

#splRep	0	1129	0
#splRep	1129	1130	3
#splRep	1130	1131	1174
print "XXXXXXXXXXXXXXXXXXXXXXXXX\n";
print "$f\n";
open(IF,"$f");
#open(IF,"test.bed");
$n1=0;$nsj=0;$nex=0;
while($a=<IF>){
$n1++;
chomp($a);
@b = split(/\t/,$a);
$end = $b[2];
$total{$end}++;
if($end<$min){$min=$end-1;}
if($end>$max){$max=$end+1;}
$hasspljun=0;
$ec = $b[9];
#print "$ec\n";

if($ec>1){
@sizes = split(/,/,$b[10]);
@starts = split(/,/,$b[11]);
for($i=1;$i<$ec;$i++){
$sjstart = $b[1] + $starts[$i-1] + $sizes[$i-1] ;
$sjend = $b[1] + $starts[$i];
#splRep	1130	1636	m150909_135227_42146_c100889372550000001823187103261633_s1_p0/5638/ccs	60	+	1130	1636	255,0,0	3	20,94,74	0,169,432
if(($sjend - $sjstart)>100){$hasspljun++;}
}
}



if($hasspljun>0){$spljun{$end}++; $nsj++;} 
elsif($hasspljun==0){$nosj{$end}++; $nex++;} 
else{print "wtf\n";exit;}
}
$n2 = $nsj + $nex;
print "(lines) $n1 nsj =    $nsj nex=$nex nsj+nex=$n2\n";
close(IF);

$min--;
$max++;
$n2=0;
print "min=$min\n";
print "max=$max\n";
foreach $end (keys %total){
#test totals!!!
#print "$t \t";
if($total{$end}){ $n2 += $total{$end}; }
}
print "$n2 $min $max\n";

#print "min=$min\n";
#print "max=$max\n";
#$n2=0;
#for($i=$min;$i<=$max;$i++){
##test totals!!!
##print "$t \t";
#if($total{$i}){ $n2 += $total{$i}; }
#}
#print "$n2 $min $max\n";
$min=1;
open(OF,">$f.dat");
open(OFSJ3,">sj3_frac.$f.wig");
print OFSJ3 "track type=wiggle_0 name=fraction_sj3_$f \nfixedStep chrom=splRep start=1 step=1\n";
$bigsjtotal = 0;
$bigtotal = 0;

for($i=$window;$i<=$max;$i++){

if( $nosj{$i}){
#ok
}
else{$nosj{$i}=0;}

if( $spljun{$i}){
#ok
}
else{$spljun{$i}=0;}

if( $total{$i}){
#ok
}
else{$total{$i}=0;}

$bigsjtotal += $spljun{$i};
$bigtotal += $total{$i};
#print "$bigtotal $total[$i] $bigsjtotal $spljun[$i] \n";
# prob need to window
# actual calc
$thistotal=0;
$thissj = 0;
for($j=0;$j<30;$j++){
$k = $i+$j + $window ;
if($total{$k}){ $thistotal +=  $total{$k};}

if($spljun{$k}){ $thissj += $spljun{$k};}

}

if($thistotal > 3){
$fracsj = $thissj/$thistotal;
}
else{$fracsj=0;}

$t = $nosj{$i} + $spljun{$i} ;
print OF "$i $spljun{$i} $t $fracsj \n";
#print "$i spl=$spljun[$i] ex=$nosj[$i] total:$total[$i] $t fsj:$fracsj \n";
print OFSJ3 "$fracsj \n";
}

close(OF);
close(OFSJ3);

system("wigToBigWig sj3_frac.$f.wig chrom.sizes sj3_frac.$f.bw ");
print "(lines) $n1 #sjtotal $bigsjtotal #bigtotal=$bigtotal \n";
#exit;
}