#! /usr/bin/perl -w
use List::Util qw[min max];

# Ashray Urs - ashray@berkeley.edu

if($ARGV[0]){
    if ($ARGV[0] eq "-calc"){
	$seqcount=0;
	calc();
    }
}

# Get information from a FASTA file
sub calc{
    if($ARGV[1]){
	$filename = $ARGV[1];
	if(-e $filename){
	    open FILE, $filename;
	    while(<FILE>){
		chomp;
		if(/>/){
		    if($sequence){
			print "There should be only one sequence in the fasta file!";
		    }
		    s/^>//;
		    $name=$_;
		    $sequence = "";
		    $len =0;
		    $seqcount++;
		}
		else{
		    $sequence .= $_;
		    $len+=length;
		}
	    }
	    dataFromSeq();
	    close FILE;
	}
    }
}

# Get and print GC content to wig file
sub dataFromSeq{
    $sequence = lc($sequence);
    @sequence = split //, $sequence;
    $filename=$ARGV[3];
    open(MYOUTFILE,">>$filename");
    print MYOUTFILE "track type=wiggle_0 name=\"GC Content\" description=\"GC Content per 5 nucleotides\"\n";
    print MYOUTFILE "fixedStep chrom=$name start=1 step=5 span=$len\n";
    for(my $i=0;$i<$len;$i+=5){
	$gc=0;
	for(my $j=$i;$j<$i+5&&$j<$len;$j++){
	    if($sequence[$j] eq "g" or $sequence[$j] eq "c"){
		$gc++;
	    }
	    $shortlen = min(5,$len-$i);
	}
	$printnum=$gc/$shortlen;
	if($ARGV[2] eq "-save"){
	    if($ARGV[3]){	    
		print MYOUTFILE "$printnum\n";
	    }
	}
	
    }
    close(MYOUTFILE);
}