#!/usr/local/bin/perl -w

#
# Generate simple gnuplot grafik (png) 
# about coverage of single/multiple contigs from ACE files
#
# Author: Bernd Senf (bsenf@fli-leibniz.de), 2008
#
# Syntax:  perl ace_contig_coverge.pl < test.ace
# 
# Output:  contig_00001.png
#          contig_00002.png
#          ....
#

use strict;

my %seq;
my %contig;

my @clen;

print STDERR "Reading ...\n";
while (<>) {
    chomp();

    if (/^RD (\S+)\s+(\d+)/) 
	{
		my ($name, $len) = ($1, $2);
		$seq{$name}{len} = $len;
		while (<>) 
		{
	    	last if (/^$/);
		}
    }

    if (/^CO (\S+)\s+(\d+)/) 
	{
		my ($name, $len) = ($1, $2);
		
		$contig{$name}{len} = $2;
		
		while (<>) {
	    	last if (/^$/);
		}

		my $qual = "";
		<>;
		while (<>) {
	    	last if (/^$/);
		}

		while (<>) {
	    	chomp();
	    	last if (/^$/);
	    	if (/^AF /) 
			{
				push(@{$contig{$name}{AF}}, $_);
	    	}
		}
    }
}

#
# sort descending
#

my $cnt=0;
my $limit=2;

foreach my $c ( sort {$contig{$b}{len} <=> $contig{$a}{len} } keys %contig )
{
	my @coverage;
	my $len=$contig{$c}{len};
	my $i;
	
	$cnt++;
	if  ( $cnt > 10 ) { print STDERR "Stop at contig number 1 ...\n"; last ; };	

	for ($i=1;$i<=$len;$i++) {
		$coverage[$i]=1;
	}
	
    foreach my $af (@{$contig{$c}{AF}}) {
		my ($name,$dir,$pos) = $af =~ /^AF (\S+)\s+(\S)\s+(\d+)/;
		my $len = $seq{$name}{len};
		for ($i=$pos; $i<$pos+$len;$i++)
		{
			$coverage[$i]++;
		}
    }
	
	my $n=0;
	my $neednewname=1;
	my $fnall="coverage_".$c.".dat";
	my @fn;
	
	open(FHALL,">$fnall") || die "Can't create $fnall:$!\n";
	print STDERR "Writing ... $fnall\n";
	print STDERR "Length  ... $contig{$c}{len}\n";

	for ($i=1 ; $i<=$len ; $i++ ) 
	{
		printf FHALL "%d %d\n",$i+1,$coverage[$i];
		if ( $coverage[$i] > $limit)
		{
			if ($neednewname) {
				my $fname="/tmp/cov_".$c."_".$n.".dat";
				open(FH,">$fname") || die "Can't create $fname:$!\n";
				$neednewname=0;
				push (@fn,$fname);
				$n++;
			}
			printf FH "%d %d\n",$i,$coverage[$i];
		}
		else 
		{
			$neednewname=1;
		}
	}
	close(FHALL);
	
	$n=$#fn;
	$i=0;
	
	if ( $n > -1 ) {
		open (PROG,"|gnuplot") || die "can't start gnuplot program ...:$!\n";
		print PROG "set terminal png medium size 1024,680\n";
		print PROG "set output \"coverage_$c.png\"\n";
		print PROG "set key off\n";
		print PROG "set xrange [1:$len]\n";
		print PROG "set yrange [1:]\n";
		print PROG "set title \"Contig $c\"\n";
		print PROG "set xlabel \"length\"\n";
		print PROG "set ylabel \"coverage\"\n";
		print PROG "plot ";

		foreach my $f (@fn)
		{
			$i++;
			if  ($i <=  $n )  { print  PROG "\"$f\" using 1:2 w i,";}
			if  ($i >   $n )  { print  PROG "\"$f\" using 1:2 w i ";}
		}
		print PROG "\n";
		print STDERR "gnuplot  => coverage_$c.png\n\n";
		
	} else {
	
		print STDERR "contig $c without coverage > $limit\n\n";
		
	}
}