#!/usr/bin/perl -w

# creates an R boxplot of omega values 
# Bruno

use strict;

$|=1;

my $TMPBOXPLOT = '_boxdata.tab';
my $RVERBOSE = 0;

my $Rparams = '';
if(!$RVERBOSE){ $Rparams = '-q 2>&1 > /dev/null' }

my ($panfile,$outfilename) = @ARGV;
if(!$outfilename){ die "# usage: $0 <tab file> <outPNGfile>\n" }

my ($maxsize,$maxoccup,$occupancy,$g,$s,$omega,%mdata) = (0,0);
open(TAB,$panfile) || die "# $0 : cannot read $panfile\n";
while(<TAB>)
{
	#<omega> pairs   n_of_seqs       file
	#0.53    1485    55      my_1787986_Bradi4g05590_1copy_aln.fna
	#print;
	next if(/^<omega>/ || !/^\d\.\d+\s+\d+/);
	my @data = split(/\t/);
	$omega = $data[0];
	$occupancy = $data[2];

	next if($omega eq 'n/a' || $occupancy < 4);

	next if($omega == 0 || $omega > 1.5);
 
	push(@{$mdata{$occupancy}},$omega);
	if(scalar(@{$mdata{$occupancy}}) > $maxsize)
	{
		$maxsize = scalar(@{$mdata{$occupancy}});
		$maxoccup = $occupancy
	}
}
close(TAB); #print "$maxsize ($maxoccup)\n"; exit;

open(TAB,">$TMPBOXPLOT");
foreach $g (4 .. 19)
{ 
	#printf TAB ("%d(%d)\t",$g,scalar(@{$mdata{$g}}));
	print TAB "$g\t"; 
} 
print TAB "\n"; 

for($s=0;$s<$maxsize;$s++)
{
	foreach $g (4 .. 19)
	{
		if(defined($mdata{$g}[$s])){ $omega = $mdata{$g}[$s] }
    #logy https://www.researchgate.net/post/Log_transformation_of_values_that_include_0_zero_for_statistical_analyses2
    #if(defined($mdata{$g}[$s])){ $omega = $mdata{$g}[$s] + 0.001 } 
		else{ $omega = 'NA'; } #$omega = $mdata{$g}->[$s] || 'NA';	
		print TAB "$omega\t";			
	}
	print TAB "\n";
}
close(TAB);

#https://sites.google.com/site/davidsstatistics/home/notched-box-plots
open(RSHELL,"|R --no-save $Rparams ") || die "# $0 : cannot call R: $!\n";
print RSHELL<<EOR;
#pdf("$outfilename",width=9);
png("$outfilename",width=800);

pancolors = c(rep('orange',14),'yellow','white')

expr = read.table(file="$TMPBOXPLOT",header=T,check.names=F); 
boxplot( expr, ylab='dN/dS',xlab='genome occupancy', notch=T, 
			varwidth=F, show.names=T, outline=F, 
			col=pancolors, cex.axis=1.5, cex.names=1.5, cex.lab=1.5 );
grid(nx=NA, ny=20,lwd=2) #grid over boxplot
par(new=TRUE)
boxreport = 
boxplot( expr, ylab='dN/dS',xlab='genome occupancy', notch=T,
         varwidth=F, show.names=T, outline=F,
         col=pancolors, cex.axis=1.5, cex.names=1.5, cex.lab=1.5 );

legend("topright", inset=.05, c("shell","soft-core","core"), 
	fill=c('orange','yellow','white'), horiz=TRUE, cex=1.5)
dev.off()

print(boxreport)

q()
EOR
close RSHELL;



